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OF POOR QUALITY 


DEVELOPMENT OF A MULTILEVEL OPTIMIZATION APPROACH TO THE DESIGN 

OF MODERN ENGINEERING SYSTEMS 


by 

Jean-Francois Marie Barthelemy 
(ABSTRACT) 

This work describes an optimization approach to the design 
of complex engineering systems. The approach assumes a 
decomposition of the original problem (design of the system) 
into smaller subproblems (design of subsystems) organized in 
a multilevel hierarchy. As one goes down the levels, the 
details of the system become more precisely defined; 
furthermore, each subproblem may control the design of lower 
level subproblems. 

A general algorithm is proposed which carries out the 
design process iteratively, starting at the top of the 
hierarchy and proceeding downward. Each subproblem is 
optimized separately for fixed controls from higher level 
subproblems. An optimum sensitivity analysis is then 
performed which determines the sensitivity of the subproblem 
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design to changes in higher level subproblem controls. The 
resulting sensitivity derivatives are used to construct 
constraints which force the controlling subproblems into 
chosing their own designs so as to improve the lower level 
subproblem designs while satisfying their own constraints. 

The applicability of the proposed algorithm is 
demonstrated by devising a four-level hierarchy to perform 
the simultaneous aerodynamic and structural design of a 
high-performance sailplane wing for maximum cross-country 
speed. The levels are devoted successively to selecting 
global performance parameters, determining the wing 
aerodynamic shape, defining the spanwise distributions of 
global structural characteristics, and, performing the 
detailed design of wing substructures. 

Finally, the concepts discussed are applied to the two- 
level minimum weight structural design of the sailplane 
wing. The numerical experiments show that discontinuities 
in the sensitivity derivatives may delay convergence, but 
that the algorithm is robust enough to overcome these 
discontinuities and produce low-weight feasible designs, 
regardless of whether the optimization is started from the 
feasible space or the infeasible one. 
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Chapter I 
INTRODUCTION 

The design of modern engineering systems is a procedure 
which often integrates different disciplines and always 
involves a large number of variables. While optimization 
techniques provide an attractive formal tool for carrying 
out the design process, a weakness of these techniques 
remains in their inability to handle truly large 
multidisciplinary design problems. To overcome that 
weakness, designers must resort to some sort of 
decomposition method. The original design problem is broken 
down into smaller subproblems which are then optimized 
separately. However, these subproblems are generally 
coupled so that an iterative scheme must be devised which 
coordinates their optimization in order that the resulting 
design is, to some extent, optimum with respect to the 
original problem. 

1 . 1 OVERVIEW OF EXISTING LITERATURE 

Several studies have been devoted to decomposition of large 
optimization problems. These studies pertain to all the 
areas of engineering as well as economics and management. 
The following overview focuses mainly on structural 
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optimization. For the sake of the discussion, two classes 
of decomposition methods are defined: formal methods and 
intuitive methods. In formal methods, the mathematical 
structure of the problem is exploited to arrive at a rigid 
decomposition scheme. Consequently, a rigorous framework 
exists within which the mathematical properties of the 
method may be assessed. In intuitive methods, understanding 
of the behavior of the physical system considered is 
apparently the prime factor directing the decomposition. 
These methods are sometimes referred to as rational methods 
and their mathematical characteristics can seldom be studied 
in great detail. An intuitive approach provides, of course, 
the only option for decomposing those problems which do not 
possess the structure for which a formal decomposition 
method exists. This division into formal and intuitive 
methods is somewhat arbitrary as a given approach may very 
well be shown to belong to both classes; however, it 
facilitates the discussion. 

1.1.1 Formal Decomposition Methods 

A very extensive body of work exists on decomposition in 
linear programming (LP). The existence of truly large 
problems in the fields of economics and operations research 
has stimulated efforts aiming at exploiting the special 
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structures of the constraint matrix. The major initial step 
in that area seems to have been the introduction of the 
Dantzig-Wolf e decomposition principle in 1960. By 1970 the 
developments in this active field were so numerous as to 
warrant the publication of a textbook by Lasdon (Ref. 1). 
In the area of structural design, problems involving 
collapse design of trusses and frames were solved 
successfully using the LP decomposition techniques (Refs. 
2-5) . 

Dynamic programming (DP) was discussed by Bellman as 
early as 1957 (Ref. 6). It is a decomposition method 
suitable for nonlinear problems. The problems must be 
serial, that is, of the form such that any change in the 
design of one subproblem affects only those subproblems that 
are "downstream" in the decision-making process. 
Furthermore, the problem objective function must be 
additive. However, DP cannot easily handle constraints that 
involve more than one subproblem. Also, it becomes very 
expensive if there are large amounts of data transmitted 
between successive subproblems. Among the advantages of the 
DP method are its capability to generate global optimum 
designs and to handle discrete variables and discontinuous 
functions. In structural design, dynamic programming seems 
to find use in optimization of "one-dimensional" structures. 
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that is, beams, one-bay multistory or one-story multibay 
frames, and transmission towers (see Ref. 7 Chap. XI or Ref. 
8 ). 

Optimization algorithms were devised for separable 
nonlinear problems using coordination techniques developed 
in the context of multilevel decision-making processes. 
Once a problem has been decomposed into smaller subproblems, 
the main task is to coordinate the design of the different 
subproblems. Essentially, the subproblems are grouped on 
the lower level of a hierarchy and an additional higher 
level subproblem is added whose role is to select the values 
of coordinating variables so as to force the other 
subproblems into choosing designs corresponding to improved 
overall problem performance. The coordinating subproblem is 
itself cast in the form of an optimization problem. As the 
evaluation of the objective function of the coordinating 
subproblem requires complete optimization of all lower level 
subproblems, this type of algorithm may turn out to be very 
expensive. A study of coordination in hierarchical systems 
is given in the 1970 monograph by Mesarovic and coauthors 
(Ref. 9). Kirsch and coworkers have used both the model 
coordination technique and the goal coordination technique 
to solve various structural design problems (Ref. 7 Chap. X 
and Refs. 10-12) 
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1.1.2 Intuitive Decomposition Methods 

The first attempts at developing intuitive decomposition 
schemes for large structural design problems were extensions 
of the fully stressed design method. In these approaches 
the structure is seen as a combination of elements 
(substructures). Given an initial design for all the 
elements, an analysis of the structure is made to determine 
inter-element forces. Then, each element is optimized 
separately on the assumption that changes in that element 
design do not change inter-element forces. Once all the 
elements have oeen designed, the structure is analyzed again 
and the process repeated until convergence is achieved. 
Giles (Ref. 13) and Sobieski (Ref. 14) performed the design 
of airplane wings under constraints on stresses and element 
stability using that approach. Kirsch and coworkers (Ref. 
15) designed frameworks using a similar approach but 
reanalyzing the structure after each substructure 
optimization, in order to account better for load 
redistribution. 

When optimizing a structure one substructure at a time, 
it is difficult to handle global constraints, that is, those 
constraints affected by variables belonging to more than one 
substructure. In Ref. 16, Sobieski and Loendorf designed 
airplane fuselages under local constraints on stresses and 
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local instabilities and global constraints on fuselage 
elastic displacements. The optimization was first carried 
out with the local constraints, as described above. If 
necessary, the resulting design was subsequently modified to 
satisfy the displacement constraints using a unit load 
method to determine the impact that changes in element 
design have on the violated displacements. Another example 
of treatment of global constraints is given by Hughes and 
coworkers (Ref. 17) who designed ship compartments for 
minimum cost under constraints on stresses. 

A generalized fully stressed design approach to large 
problems is certainly economically appealing. However, it 
presents two difficulties. As pointed out in Ref. 14 
"...minimization of the individual component masses does not 
guarantee minimization of the total mass . This situation 
is caused by the inability to control the load path on the 
assembled structure level... ." As mentioned earlier also, 
it makes it difficult to handle global constraints. In 
Refs. 18 and 19 Schmit and coworkers used a two- level 
approach to design trusses and aircraft wings. At the 
global level, the distributions of stiffnesses, the global 
level variables, were chosen so as to minimize the total 
structural weight, while satisfying global displacement 
constraints and also some local constraints on stresses and 
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structural element buckling. For known stiffnesses, the end 
forces on the various structural elements were calculated. 
At the local level, these structural elements were optimized 
separately with respect to their detail design, the local 
level variables, so that the changes in element stiffnesses 
were minimized, while the local constraints were still 
satisfied. The process was repeated until convergence was 
achieved. In this decomposition, the introduction of the 
global level problem was a key factor in overcoming both 
difficulties attributed to the generalized fully stressed 
design approach. By placing the minimization of the weight 
at the global level, the opportunity was kept to trade 
structural mass between the elements in order to improve the 
load paths while reducing the total weight. Also, the 
provision was left to explicitly handle global constraints. 

1.1.3 General Approach to Engineering System Design 
The intuitive decomposition methods discussed to this point 
are very specific and cannot easily be extended to different 
classes of problems. In particular, they are not suited for 
multidisciplinary problems. In Ref. 20, Sobieski proposed a 
general approach to the design of large engineering systems, 
and in Ref. 21, he proceeded to demonstrate the proposed 
concepts by applying them to the two- level minimum weight 
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design of a portal framework under local constraints on 
stresses and local buckling and global constraints on 
displacement. At the global level, the weight of the 
framework is minimized with respect to the cross-sectional 
area and bending moment of inertia of the different beams, 
while satisfying constraints on displacements. Once a 
global design is obtained, the framework is analyzed to 
obtain the end loads on the beams. At the local level, each 
beam is designed separately. Its detailed dimensions are 
proportioned so as achieve the values of cross-sectional 
area and bending moment of inertia chosen at the global 
level while satisfying constraints on stresses and local 
buckling. This is done by minimizing a measure of the 
violation of the subproblem constraints. An optimum 
sensitivity analysis (Ref. 22) is then performed in order to 
obtain a linear approximation to the constraint violation of 
each local level subproblem in terms of the global level 
variables. The whole optimization is repeated, but, at the 
global level, it is required that the (linearized) local 
level constraint violations be reduced. This iterative 
process is pursued until all the local and global level 
constraints are satisfied and the framework weight cannot be 
reduced further. Although the proposed scheme is applied to 
a structural synthesis problem, it is applicable to any type 
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of engineering design. Furthermore, it may be used with 
decompositions involving more than two levels. Some of the 
key features of the approach are discussed below. 

The method postulates a multilevel decomposition of the 
problem, that is, a decomposition where the different 
subproblems are grouped on the levels of a hierarchy. The 
subproblem at the top level is concerned with the 
optimization of the overall system objective in terms of 
design variables describing that system in very global 
terms. As one proceeds down the hierarchy, the subproblems 
are dealing with subsystems described in deeper and deeper 
levels of detail. The choice of the specifics of the 
decomposition for a given application still remains largely 
a matter of physical insight and convenience. 

Below the highest level, each subproblem is concerned 
with finding a set of design variables which satisfy its own 
constraints for fixed values of higher level variables. It 
turns out that this amounts to finding a design vector which 
satisfies a set of equality constraints and a set of 
inequality constraints. Such a mathematical problem may not 
have a solution. It may be transformed into a mathematical 
programming problem which maximizes a measure of "how well" 
the design satisfies the sets of constraints. The 
subproblem formulation of Ref. 21 uses the equality 
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constraints to eliminate dependent variables and forms a 
penalty function with the inequality constraints. That 
penalty function is a measure of constraint violation which 
is then minimized with respect to the independent variables 
to find that design which best meets the constraints. 

Once designs are obtained for all the subproblems on a 
given level, sensitivity analyses are performed to determine 
the effect changes in higher level variables have on lower 
level constraint violations. The resulting sensitivity 
derivatives are used to construct linear approximations to 
the lower subproblem constraint violations in terms of the 
higher level design variables. These approximations are 
added as constraints for the next optimization of the higher 
level subproblems. This enables the optimization of the 
higher level subproblems to be conducted so as to improve 
the design of the lower level subproblems which have 
violated constraints or so as to reduce the margin of 
satisfaction for those subproblems for which all constraints 
are amply satisfied. In essence, this makes the trade-offs 
between the different subproblems visible and allows them to 
be resolved in a rational fashion. 

The task of the highest level subproblem remains the 
optimization of the system performance index with respect to 
its own global variables which are required not only to 
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satisfy its own constraints but also to reduce the 
constraint violation for the lower level subproblems. 


1.2 OUTLINE 

The purpose of this work has been to investigate the 
applicability of Sobieski's method to the design of complex, 
multidisciplinary engineering systems. It was elected to 
focus the effort on a specific example: the simultaneous 
aerodynamic and structural design of a high-performance 
sailplane wing for maximum cross-country speed. It was 
determined that this design problem could be carried out by 
means of a four-level decomposition. An algorithm was 
developed to solve the problem and it was tested on the two 
lower levels of the decomposition. 

In Chapter II, we discuss the optimum sensitivity 
analysis technique. We show that the stationary conditions 
satisfied at a local constrained optimum may be used to 
determine the sensitivity of optimum objective functions and 
design variables to changes in those parameters that were 
kept fixed during the optimization. 

Chapter III is devoted to a detailed description of the 
proposed algorithm. It is constructed along the same lines 
as that described in Ref. 20. The main difference resides 
with the specifics of the mathematical formulation for the 
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individual subproblems. The present algorithm treats 
equality constraints in a general numerical fashion. 
Furthermore, it is designed as a constrained programming 
problem rather than as a penalty function minimization 
problem. Particular attention is paid to the potential 
difficulties associated with decompositions that involve 
more than two levels. 

In Chapter IV, the decomposition proposed for the example 
problem is detailed. It is a four-level decomposition. The 
highest level subproblem entails selection of global 
performance parameters which include weight and aerodynamic 
characteristics. Then follows an aerodynamic subproblem 
where the shape of the wing is defined so as to obtain the 
aerodynamic characteristics specified in the performance 
subproblem and control the initiation of stall. On the 
third level of the decomposition is a global structural 
subproblem where the wing spanwise distributions of weight 
and stiffnesses are selected for fixed total weight and so 
as to limit the wing bending response and torsional 
divergence speed. Finally, at the lowest level, the wing is 
modeled by a series of spanwise elements each of which is 
designed separately for fixed weight and stiffnesses under 
stress and local buckling constraints. 
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In Chapter V, the algorithm introduced in Chap. Ill is 
applied to the two lower levels of the hierarchy described 
in Chap. IV. The problem solved is the minimum weight 
design of a straight composite wing. Results are given for 
designs started from both the feasible space and the 
infeasible space. 

Finally, Chapter VI is devoted to synthesizing the 
experience gained with this work and making recommendations 
for continued investigation. 
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SUMMARY OF PERTINENT RESULTS FROM OPTIMUM 
SENSITIVITY ANALYSIS 

Optimum sensitivity analysis is a technique which permits 
one to investigate the sensitivity of the solution of an 
optimization problem to variations in the parameters of the 
problem. It yields derivatives of the optimum values of the 
design variables and objective function with respect to the 
parameters, derivatives which may then be used to perform 
trade-off analyses. These derivatives are called 
sensitivity derivatives to distinguish them from the 
derivatives of objective function and constraints with 
respect to the design variables which are termed behavior 
derivatives. This technique was recently introduced by 
Sobieski and coauthors (Ref. 22) as a tool for structural 
synthesis. It has been specialized by Schmit and Chang 
(Ref. 23) to a formulation of the problem of structural 
sizing for mimimum weight based on approximation concepts 
and dual methods, and has been used by Haftka (Ref. 24) in 
the optimization of damage tolerant structures. The 
elements of the technique required for the remainder of this 
work are discussed in this chapter; the reader interested in 
a more general discussion should consult Ref. 22. 
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2.1 STATEMENT OF THE PROBLEM 

Assume that we start from the following nonlinear 
mathematical programming problem 


min F(X,P) 

X 

so that g(X,P) 5 0 (2.1) 


where vector X contains the n design variables of the 

problem, while vector P contains k parameters that are kept 

fixed during the optimization process. The constraint g is 

vector-valued, in general; it is assumed here to have m 

components. The optimum solution is functionally dependent 

* 

on the parameters, hence, using ( ) to denote optimum 

quantities 

* * 

X = X (P) 

F* = F*(P) (2.2) 


At the constrained minimum, the following stationary (Kuhn- 
Tucker) conditions hold true. 


F 

g 


* 

X 

a 




X 


★ 


0 


* 

X > 0 


(2.3) 
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In these equations, subscript x indicates partial 

derivatives with respect to the design variables; 

superscript a refers to the constraints exactly satisfied 

(active) at the local minimum (there are m such constraints 

'a 

and g is a subvector of g); finally, X is the vector of 

Lagrange multipliers (dual variables) . The object of 

optimum sensitivity analysis is to find the rates of change 
* * 

of F and X with respect to P. We define these rates of 
changes as follows 


X' 

* 

dX 

(matrix. 

nxk components) 



dP 

* 




F' 

_ dF 

(vector. 

k components) 

(2.4) 


dP 




If 

F' and 

X' are 

known, we may construct 

linear 


approximations to the solution of Problem (2.1), 

F*(P+AP) = F*(P)+F'^AP 

X*(P+AP) = X*(P)+X''^AP (2.5) 

The two following sections are devoted to calculating F' and 
X' . 


2.2 SENSITIVITY DERIVATIVES OF THE OBJECTIVE FUNCTION 
The optimum sensitivity derivative of the objective function 
is obtained by the chain rule of differentiation as 


F' 


( 2 . 6 ) 
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where subscript p indicates partial derivatives with respect 
to vector P. Now, we can show that F' may be found without 
actually calculating X’ . We require that the optimum point 
for P+AP remain at the intersection of the constraints 
determining the optimum point for P. Hence, the second of 
Eqs. 2.3 is stationary with respect to P, 


U/d\ d ^ 

ap(g ) = 9p X = 0 


(2.7) 


Combining Eq. 2.7 transposed and postmultiplied by X and 


the first of Eqs. 2.3 premultiplied by X 


.T 


we obtain 


. T * a *T * 

X'-^F = g® -^X 

X ^p 


( 2 . 8 ) 


or 


. * a*T * 

F' = Fp+g“ ^ (2.9) 

* 

The Lagrange multipliers X may be available as by-products 
of the optimization scheme; otherwise, they may be obtained 
by QR decomposition or from the more classical least-squares 
solution to the first of Eqs. 2.3, 


* 

X = 
Then, 


, a* a*T,-l a*_ 


Eq. 2.9 yields F' 


directly . 


(2.10) 

An alternate solution 


exists because Eq. 2.7 is the only condition that X' must 
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meet in order for Eq. 2.8 to be satisfied. Therefore, if 
for some matrix Y (nxk) 




Y 


0 


( 2 . 11 ) 


then 




* 


T 

Y F 


* 

X 


and 


F’ = F*+yV (2.12) 

Equation 2.11 represents k sets of m linear equations in n 
unknowns, with m S n. We may therefore give arbitrary 
values to n-m components of each column of Y and deduce the 

d. 

remaining m components from Eq. 2.11. In general, Y 
differs from the true sensitivity derivatives of the 
variables X' , unless n = m . Assume that, in the space of 
the design variables X, an infinitesimal perturbation of 
parameter P results in a move from the initial point in the 
direction specified by Y. Then, the resulting design is 
still at the intersection of the originally active 
constraints because condition 2.7 is enforced; however, it 
may not be optimum because the first of Eqs. 2.3 may not be 
satisfied any more. 
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2.3 SENSITIVITY DERIVATIVES OF THE DESIGN VARIABLES 
The first two of Eqs. 2.3 must remain satisfied at P+AP. 
Defining the sensitivity derivative of the Lagrange 
multipliers 


X- = ^ 

dP 


(2.13) 


we have 


d , a* . 

dp'9 ) 


g® +g^ X' 


= 0 


and 


d * a*T * 

^(F ) 

dP' X ' 


* * . *T a* *T a* a*T 

= F +F X'+X g +X -^g® X'+g 
xp XX ^xp ^xx ^x 


= 0 


This yields the sensitivity equations, k systems of n+m, 

c 

linear equations in n+m unknowns. In matrix form, we have 

a 


* *T a * 
F....+X -^g^ 


XX 


XX 


a*T 




\ 

/ 


X' 

1 > 

< 

X' 

< 


1 

J 

\ 


* *T a* 

F +X g 
xp ^ xp 


(2.14) 


Once Eq . 2.14 is solved for X' , Eq. 2.6 gives F' directly. 
Reference 22 discusses in detail the solvability of the 
sensitivity equations. 
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The approach of Sec. 2.2 yields the exact value of F' at 
considerably lower cost than that of Sec. 2.3 since it 
involves only first-order partial derivatives of F and 
with respect to the design variables and parameters of the 
problem. However, it does not usually permit one to obtain 
X' (unless n = m ) . The choice between the two types of 
analysis is very much dependent on the application 
considered and the nature of the optimization problem at 
hand. It must be guided by a comparison between the cost of 
calculating second-order partial derivatives of F and g on 
the one hand, and the benefit of knowing the sensitivity 
derivatives of the variables on the other hand. 

It must be noted that higher order sensitivity 

* * * 

derivatives may be obtained for X , F and X by setting 
higher order total derivatives of the two first stationary 
conditions (Eqs. 2.3) to zero. In particular. Ref. 25 shows 
how little additional calculation is required to find F”, 
once X' has been obtained. 

2.4 ADDITIONAL COMMENTS 

After solving a mathematical programming problem and 
before performing sensitivity analysis, it is necessary to 
find the set of the active constraints. These constraints 
must satisfy Eqs. 2.3; furthermore, their gradients must be 
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linearly independent, if one wants to obtain the Lagrange 
multipliers from Eq. 2.10. The algorithm described below 
was used with success. 

i) Select an initial set of nearly satisfied 

constraints by retaining those constraints whose 

value is above a small negative tolerance. 

' a a 

ii) Identify S , the subset of S which contains the 

maximum number of constraints with independent 
gradients. To achieve that, the constraint 

gradients are taken as the columns of a rectangular 
matrix and a Gaussian elimination with row and 
column interchanges is performed to identify the 
independent columns. 

iii) Use Eq. 2.10 to find the Lagrange multipliers 

' a 

corresponding to set S . Note that the fact that 

' a 

the constraints of set S have independent 
gradients is a guarantee that Eq. 2.10 has a 
solution. 

iv) Find the smallest Lagrange multiplier 

v) If is negative (less than a small positive 

tolerance), eliminate constraint i from set and 
restart the process in ii). If is positive, 

* a 

terminate the process since S is the desired set. 
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Note that steps ii) and iii) would not be required if QR 
decomposition were used. In step v) , only the constraints 
with strictly positive Lagrange multipliers are kept, even 
though active constraints may have Lagrange multipliers that 
are zero. This is because constraints with zero Lagrange 
multipliers are redundant; hence, their gradients are 
linearly dependent on other constraint gradients. 

In this work, the sensitivity derivatives are used to 
extend optimization information over a range of variation 
for the parameters of the design problems considered. This 
is achieved by means of the linear extrapolations given in 
Eq. 2.5. The implicit relationships 2.2 are expected to be 
nonlinear, in general. Hence, the extrapolations being 
first-order Taylor series expansions, they are valid only 
for limited AP. Furthermore, the calculation of the 
sensitivity derivatives is based on the assumption that the 
set of constraints active at the initial optimum will remain 
unchanged as the parameter P is changed to P+AP. 
Perturbations of the parameters cause moves of the 
constraint boundaries with respect to each other and, 
therefore, may very well result in changes in the active 
constraint set. These changes then imply reduced accuracy 
of the extrapolated optima. The question arises of whether 
it is possible to estimate the magnitude of the parameter 
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perturbations which cause changes in the active constraint 
set. It was proposed in Ref. 22 that the perturbation of 
parameter p (element of vector P) which results in 
constraint i changing from active status to inactive status 
may be obtained to the first order by stating that is 

reduced to zero when constraint i leaves the active set, 

^*(P+Ap^) = X*(p)+X^(p)Ap^ = 0 

or 

APi = -X*(p)/X^(p) (2.15) 

Reference 22 likewise suggested that the- perturbation of 

parameter p causing constraint i to change from inactive 

status to active status could be found as that value Ap^ for 
* 

which gj^(p+Ap^) becomes zero. 

g^(p+Ap^) = gj^(p) + [g*^(p)X' (p) ]Ap^ = 0 
then 

APi = -g*(p)/Ig*^(p)X' (p) ] (2.16) 

It was observed in Ref. 26 that Eqs . 2.15 and 2.16 could be 
unreliable. These estimates were used to determine with a 
specific design problem example what parameter perturbations 
cause changes in the status of the different constraints 
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involved, attempting , so to speak, to construct an history 
of the active constraint set as p varies. This may not be 
an appropriate approach. Indeed, Eqs . 2.15 and 2.16 may be 

used to estimate a Ap^ for each of the problem constraints, 
indicating when the status- of a particular constraint is 
supposed to change. However, only the smallest (in absolute 
value) negative and positive of these Ap^ (Ap_ and Ap^ ) 
should be taken into consideration. The extrapolations 
based on sensitivity derivatives obtained at p should be 
adequate for p+Ap_ S P+Ap < p-+Ap^ . But any prediction 
beyond that range (in particular prediction of change in the 
active constraint set) should not be relied upon since it 
would have been obtained with derivatives corresponding to 
an inappropriate active constraint set. In essence, as 
suggested in Ref. 23, the Ap^ of Eqs. 2.15 and 2.16 must be 
regarded as bounds on Ap. Barring excessive nonlinearity, 
the range of Ap defined by the intersection of these bounds 
is that within which the extrapolations based on sensitivity 
derivatives calculated at p should be reliable. In the 
absence of changes in the active constraint set, the 
accuracy of Eq. 2.5 is limited only by the inherent 

nonlinearity of the specific design problem considered. 

Experiments conducted with various structural design 
problems showed (Ref. 22) that linear-extrapolation-based 
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variables fell within a few percent of actual optimization 
result, for parameter variations of up to twenty percent; 
the predictions often were better for the objective function 
than for the design variables. 
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Chapter III 

PRESENTATION OF THE ALGORITHM 

This chapter is devoted to discussing an algorithm proposed 
for the design of complex engineering systems. The method 
is based on multilevel decomposition of the corresponding 
design problems into subproblems which are solved 
separately. The algorithm is discussed in general terms. 
This is done so as to avoid focusing on a particular 
problem. The following chapter will be devoted to 
investigating the applicability of the approach described 
here to a specific problem. The abstract concepts 
introduced here can be more easily put in a physical context 
if this chapter and the following one are read 
simultaneously. The presentation begins with a definition 
of the type of fhultilevel decomposition considered and of 
the components of that decomposition. Then, the discussion 
centers on the formulation and solution of the subproblems. 
Finally, the organization of the overall optimization is 
discussed. The developments follow those by Sobieski in 
Ref. 20. 
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3 . 1 DEFINITIONS 

The starting point of this discussion is the task of 
designing a system by finding a set of variables which 
optimize a performance index while satisfying a number of 
constraints. This will be referred to as the 

overall problem . Assume that, for convenience, the overall 
problem is decomposed into a number of subproblems , each of 
which corresponds to the design of a subsystem, and has its 
own sets of variables and constraints. The subproblems are 
organized in a hierarchy as shown on Fig. 1. The subproblem 
at the top of the hierarchy is the design of the whole 
system, in very general terms. The next level deals with 
major subsystems, the third one with subsystems of the major 
subsystems, and so on. As one goes down the levels, the 
details of the system become more and more precisely 
defined. The subproblems are identified by two numbers ij 
which define (i) the level of the subproblem, and (j) its 
position on that level. The level numbers increase 

downward, and the subproblem numbers increase, say, from 
left to right. The hierarchy is assumed to have 1 levels; 
there are n^ subproblems on level i . 

The whole design process is conducted iteratively, 
starting at the highest level and going down. Each 
subsystem is designed separately. The object of the 
algorithm is to coordinate the design of the different 
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Legend: 



subproblem ij 


response 


Figure 1: 


Example of four-level hierarchy. 
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subsystems so as to obtain optimum system performance by 
resolving the trade-offs between the subsystems in a 
rational fashion. This is achieved by manipulating two 
types of interactions between the subproblems: the controls 
and the responses. A control is a quantity selected in a 
subproblem to influence another subproblem in the hierarchy. 
A response is a quantity that reflects the effect the design 
of one subproblem has on that of another subproblem. In 
this work, we assume that the controls flow exclusively 
downward in the hierarchy, while the responses flow 
exclusively upward. This excludes situations where a 
subproblem controls another subproblem on the same level 
( lateral control ) or on a higher level ( reverse control ).^ 
As represented symbolically on Fig. 2, each subproblem is 
assumed to receive controls from higher level subproblems 
through vector (control inputs in subproblem and 
to receive responses from lower level subproblems through 
vector (response inputs in subproblem ij) . Also, each 
subproblem is assumed to control lower level subproblems 
with vector (control outputs from subproblem ^), and 


In Ref. 20, the controls are named interaction quantities. 
The interested reader should refer to that source for 
indications as to how to handle reverse and lateral 
controls . 


1 
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c) Highest level subproblem. 


i 

RO. . 

11 

1 

Cl . . 

11 

1 


Subproblem ij 

variables X. . 

11 

initial design 

x°. 

11 

i 

RI . . 
11 

l 

CO. . 
11 

1 

r 


a) Lowest level subproblem. 



Intermediate level subproblem. 



Figure 2: 


Symbolic representation of the different types of 
subproblems . 
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to respond to higher level subproblems with vector 

(response outputs from subproblem ij^). This discussion does 

not make any assumption on the specific form of the 

controls, except that the control outputs from a subproblem 

depend on the variables of that subproblem; they may be the 

variables themselves or functions of these variables (CO. . = 

C0^j(X^^)). Likewise, the controls imposed on subproblem ij 

(Cl^j) are not specified , except that they are assumed 

known once all the control outputs from higher level 

subproblems o.<±) are known. The response outputs 

from a subproblem (RO^^) will be defined in this chapter. 

Furthermore, it will be shown that the response inputs into 

a subproblem (RI^^. ) are entirely known once the response 

outputs from the lower level subproblems (RO „, o>i) are 

ap 

known . 

For each iteration of the overall system, each subproblem 
is concerned with finding a design which satisfies its own 
constraints for the controls specified by the higher level 
subproblems (Cl^j). In addition, the subproblem manipulates 
its own controls (00^^ ) so as to improve the responses from 
the lower level subproblems (RI^^). Mathematically, this 
implies finding a design vector that satisfies a set of 
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equality constraints and a set of inequality constraints. 
This problem may be solved by transforming it into a 
nonlinear mathematical programming problem where a measure 
of constraint violation is minimized with respect to the 
design vector. The subproblem's own response (RO^^) is then 
constituted of the minimum measure of constraint violation 
as well as derivatives indicating the sensitivity of that 
minimum measure to changes in the controls imposed on the 
subproblem. 

The subproblem on the highest level of the hierarchy 
drives the whole design process. It selects its own design 
variables so as to optimize the overall problem performance 
index and satisfy its own constraints. It adjusts its 
controls (CO^^) to improve the responses from the lower 
levels (RI^^), that is, reduce the constraint violation in 
the lower level subproblems. 


3.2 LOWEST LEVEL SUBPROBLEM 


Figure 2 a) describes symbolically subproblem 1 j , a 
subproblem on the lowest level of the hierarchy. The vector 
of design variables is , the subproblem design is 
influenced by the controls collected in vector CI^^ and the 
subproblem response is vector subproblem 


statement is as follows. 
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(3.1. a) 



(3.1.b) 



(3. l.c) 


The functions and are vector-valued. The subproblem 


design is subject to a set of inequality constraints (Eq. 
3.1. a). These are the constraints set on the corresponding 
subsystem in the statement of the overall non-decomposed 
problem. The design is also subject to a set of equality 
constraints (Eq. 3.1.b). Some of these may come from the 
formulation of the overall non-decomposed problem; some are 
imposed on the subproblem during decomposition. The 
multilevel decomposition consists of different models of the 
system being designed, the models becoming more refined as 
one progresses down the hierarchy. The equality constraints 
arise from the necessity of requiring consistency between 
the models at different levels. Also, some upper and lower 
bounds may be set on the variables (Eq. 3. l.c). As stated 
in the previous section, the subproblem is assumed to be 
controlled by each subproblem on the higher levels. As a 
result, the vector of controls functionally 
dependent on the optimum designs of all the subproblems 
above the lowest level. 
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Cli. = 



(3.2) 

Problem 3.1^ may have from 

zero to 

an 

infinity of 

solutions, depending on the controls 

We 

transform it 

into a nonlinear mathematical 

programming 

problem by 

introducing two new variables 

and which 

are used to 

relax inequality and equality 

constraints . 

These new 

variables are combined into a composite function which is a 

global measure of constraint 

violation 

and which is 

minimized to identify the "best" 

design . 

The 

transformed 

problem is 




find j ' so that 

-E , . 

e is minimum 



(3 .3 .a) 




(3.3.b) 

-"ij ^ ^ "ij 



(3.3.C) 




(3 .3 .d) 

1 

lA 

O 



(3.3.e) 


^ In the remainder of this dissertation, references are made 
to "Problem 3.j", meaning the mathematical problem 
described by Egs. 3.j. 
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a) Initial problem has feasible solutions (e* k 0). 


re 3: Illusrration of the principle of inequality 

constraint relaxation. 
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where b is a positive constant. Each equality constraint is 

replaced by two inequality constraints (Eq. 3.3.c). The 

distance between these constraints (2 ti^^) is forced to 

remain positive because of Eq. 3.3.e and to decrease because 

of the second term in the objective function. Each 

inequality constraint is relaxed by the addition of variable 

E, .. The mechanism of that relaxation is described in Fig. 
Ij 

3, for a simple two-dimensional problem with one equality 

and three inequality constraints, assuming that the equality 

constraint is identically satisfied. The first term in the 

objective function (Eq. 3. 3. a) is defined so as to increase 

Ej^j. A high positive value for e^^ is quite acceptable as 

it means that the selected design point is deep inside the 

feasible domain, away from the constraint boundaries, and 

that decreases the chances of obtaining infeasible designs 

when the controls are perturbed. If the solution of Problem 

3.3 is so that e, . SO and ti, . = 0, then it solves Problem 

Ij l3 

3.1. If, on the other hand, we have e^^ < 0 or » 0, 

then the resulting design is infeasible with respect to 

Problem 3.1. 

Constant b is added to the definition of the objective 
function to control the importance attributed to violating 


the different constraints. A high value for b results in a 

heavy penalty paid for violating equality constraints; 

conversely, a low value for b makes the penalty for 

violating inequality constraints more important, in relative 

terms. In addition to adjusting b, it may prove appropriate 

to scale variables e, . and ti, . in order to make all 

Ij '1: 

components of the constraint gradients of the same order of 
magnitude . 

Any nonlinear programming algorithm may be used to solve 

★ 

Problem 3.3. This yields ^ , as well as a measure of 

* 

inequality constraint violation and a measure of 

* 

equality constraint violation . These are functionally 

dependent on the controls and, therefore, as indicated 

by Eq. 3.2, on the optimum design variables of subproblems 

on higher levels. Optimum sensitivity analysis may be used 

to study the effect changes in higher level subproblem 

variables have on the constraint violations for subproblem 

* 

Ij . Considering if higher level subproblem 

designs are perturbed, we have 


e! JCI, . (X* +AXt ,x! T +AX. ^ 

I3 I3' 11 11' 1-1, n^_^ 1-1, n^_^ 


)] 


1-1 n 


■ , , 


(3. 4 . a) 


where 


* 

^Ijop 
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,X 


1-1 , n 


)] 


1-1 


+ 


3e, . 9CI, . 

Ij l2 

3CI, . 3X 
Ij a3 


(3.4.b) 


and 


0 < T 


IjaP 


< 1 


(3.4.C) 


1-1 n 

S t ^ 

Equation 3.4 is a first order Taylor series expansion of the 
type discussed in Chap. 2. It indicates the sensitivity of 
the design of subproblem Ij to changes in higher level 
subproblem designs. It is written as a superposition of 
responses to each of these subproblems. Equation 3.4.b 
gives the response to subproblem o&. Several subproblems 
may control subproblem 1 j . Hence, each one is given the 


task of reducing a fraction of The first term in Eq. 

3.4.b determines the fraction of constraint violation in 

subproblem Ij that must be reduced by appropriate choice of 

the controls of subproblem a&. The coefficient Taay be 

defined as the coefficient of participation of subproblem cP 

in the design of subproblem Ij . It must be chosen on the 

basis of experience with 

* 

Because of Eq. 3.4.d, 


the type of system considered, 
is entirely divided among all 


39 


PhOZ 

OF POOR QUALITY 


subproblems controlling subproblem 1 j . The second term in 
Eq. 3.4.b indicates how changes in the design of subproblem 
a& will affect the design of subproblem Ij . The same 

•k 


linearization may be conducted for 
comparable information for the measure 
equality constraints. 

The vector of response outputs from 
defined as 


providing 
of violation of 

subproblem Ij is 



< 


* ^^^Ij * 

i 


3n 


Ij 


3CI 


Ij. 


3CI 


IjJ 


(3.5) 


This information will be used to construct the response 
inputs into higher level subproblems «<1). 


3.3 INTERMEDIATE LEVEL SUBPROBLEM 

Figure 2 b) depicts symbolically a subproblem on any 
intermediate level of the hierarchy. The design variables 
are arranged in vector . At the beginning of the present 
iteration, the initial value of is ^ , the optimum 

design obtained during the previous iteration. The design 


of the subproblem is influenced by the controls CI^^ and the 
responses The subproblem influences the lower level 
subproblems with its own controls ^ ; its response to the 
higher level subproblem controls is RO . .. The controls Cl. . 

ij ij 
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are issued £rom all tlie higher level subproblems. 

the following functional relationship 

We have 

Cl . . = Cl . . (X* , . . . ,X* T ) 

ij i:' 11' 1-1, n^_^' 

(3.6) 


The definition of will be discussed at the end of this 


section. The statement of the subproblem is 
find X. . so that 


g. . (X. . ,CI . . ) 1 0 

^ij ' 1 j ' 1 j ' 

(3. 7. a) 

e. . (X. . ,CI . . ) =0 
ij ij 13 

(3.7.b) 

-t* . .(AX. .) = -£* . .(X. .-X°.) < 0 
yvij ' ij ' yvij ij ij 

(3.7.C) 

T 1 * . . (AX. . ) = Ti* . . (X. .-X° . ) = 0 
vvij ' X2 yvij ' 13 13 ' 

(3.7.d) 

V = i+l,...,l; V = l,.,.,n^ 

X^- . 5 X . . 5 x'^ . 

13 13 13 

(3.7.e) 

(1-6)X° . S X. . ^ (1+6)X° . 

' 13 13 13 

(3.7.f) 


This statement contains three types of constraints which 
were not included in the formulation of the lowest level 
subproblem. Constraints 3.7.c and d are linear 
extrapolations of the type of Eq. 3.4.b, their coefficients 
being available from the responses j • They require 
satisfaction of equality and inequality constraints in lower 
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level subproblems. Equation 3.4.b gives V = 1; a 

more general expression is given at the end of this section. 

The increment corresponds to the change of with 

respect to the previous optimum design, or, the current 

initial design, hence AX. . = X. .-X?.. Constraints 3.7.f 

ij ij 13 

are called move limits; they limit the changes in the 
variables X^^ to a fixed fraction 6 of the initial values of 
these variables. This is done to preserve the validity of 
the sensitivity-analysis-based linear extrapolations of Eqs . 
3.7.C and d (see Sec. 2.4). 

The solution proposed for this problem is very similar to 
that used for Problem 3.1, namely: 


find X. .,E. .,Ti. . so that 
13' 13' 13 


-e 


11^1- • • 

e -^+bTi^^ is minimum 


g. . (X. . ,CI . . )+e . . ^ 0 

^13' 13 11^ 13 


-r) . . ^ e. . (X. . ,CI . . ) ^ T) . 
13 13' 13' 13^ '1; 


-E . . (X. .-X . )+E . . < 0 

yvi3' 13 13^ 13 


-n . . s n . . (X. .-X . ) ^ T) . . 
'13 W13' 13 13' '13 


(3. 8. a) 
(3.8.b) 
(3.8.C) 
(3 .8.d) 
(3.8.e) 


V = i + l,...,l; V = l,...,n. 


-Ti . . ^ 0 
11 

x^ . < X. . s x':^ . 
11 11 11 


(3.8.f) 


(3.8.g) 
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(3.8.h) 


The comments relative to Problem 3.3 apply to Problem 3.8 

as well. Solution of Problem 3.8 using any nonlinear 

programming algorithm yields an optimum design an 

★ 

optimum measure of inequality constraint violation 

and an optimum measure of equality constraint violation 
* 


Ti^ j . These are now functionally dependent on the controls 
Cl . . , the responses RI . . , and the initial design X° . . 

ij ^ ij' ^ ij 

Taking, for example, the measure of violation of inequality 
constraints, we have 


* * * 

t. . = E . . [Cl . . (X 
ij 13 13 


11 ' 


* 

'^i-l,n^_^' -ij '”ij 


),RI, ,,X° ] 


(3.9) 


Again, we may resort to optimum sensitivity analysis to 
investigate the effect that changes in higher level 
subproblem designs have on the design of subproblem i j . 
However, we must also account for the changes in the lower 
level subproblem responses and in the initial design for 
subproblem i j . 


: . . |CI . . (Xt T +AX, ,X° .+AX° . ] 
ij ij'll 11 ' ij ij 


i-1 n 



^ijaB 




(3. 10. a) 


where 
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E . . „(AX „) = . [e .. [Cl .. (X^ ^ / X . ] 

lja3 ija& I IJ IJ ' 11 ' IJ 


* 

3e . . 
IJ 


3RI . . 
iJ 


* 

3e . . 
IJ 


ARI . . + 
IJ 


and 


3X 


-AX°^ 


iJ 


3e . . 3CI . . 
11 11 


-AX 


3CI . . 3X „ 
11 a& 


o3 


(3. lO.b) 


0 5 4'. . „ ^ 1 

1 ja& 


(3. lO.c) 


i-1 n 

Z / f = 1 


(3. lO.d) 


This equation corresponds to that developed for a lowest 

level subproblem (Eg. 3.4). The only difference resides in 

the addition of terms accounting for perturbations in 

response inputs (ARI^^) and in initial design (AX?^). A 

* 

similar expression may be developed for ^ . 

The vector of response outputs from subproblem ij is 
defined as 


RO. . 
11 




f * * * 

dz . . 8e . . 3e . . 

11 11 11 

* * 

3CI . . 3RI . . 3X ^ 


11 


11 11 


9ti • • 

7^ 

3t1 . . 

9ti . . 

11 

11 

IJ 

> 

3CI . / 

3RI . . 

' 3X° . 

11 

11 

IJ J 


1 < i < 1 


(3.11) 


The vector of response inputs to subproblem ij contains 
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the information necessary to construct the linear 
approximations of Eqs . 3.7.c and d. 


T I * T * T’ * 'T 

^^ij ^ I'^i + l^lij^ ' • • • ' ^'^In^ij ^ '®^i + l,lij^ '• 


; T| , . . } 

In^ij 




yvij 


dc 

= . . (e +_Ji^ ARI - 

yvij \ yv yv 






3RI 


yv 


3X 


yv 


l<i<l; i<y (3.12) 

* 

{n . .} is defined similarly. There are two subvectors 

yvij 

for each lower level subproblem yv. Their components are 
obtained recursively from Eq. 3.10.b, redefining ij as yv, 
and aB as i j . To calculate the response inputs into 

subproblem ij, we need to know, for each lower level 
subproblem yv, the response outputs (RO^^), the changes in 
response inputs (ARI^^), the changes in initial design 
(AX°^), and the specific form of the relationship between 
the control inputs in subproblem yv and the variables of 
subproblem ij (CI^^(X^^)). The contributions of a 

subproblem on the lowest level ( ^^ivi j ^ ' ^’’ivi j ^ ^ 
involve terms in ARI^^ and 4X°^, as there are no response 
inputs nor move limits for a lowest level subproblem (see 
Problem 3.3). 


3e* 3CI I 
yv __H 

3CI 3X. . I 
yv ij J 
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3.4 HIGHEST LEVEL SUBPROBLEM 

The highest level subproblem is depicted in Fig. 2 c). Its 
task is to select the variables which optimize the 

overall problem objective while satisfying a set of 
inequality constraints. In addition, the controls must 

be manipulated according to RI 22 3.12) so as to improve 

the lower level subproblem constraint violations. Assuming 
that the overall problem objective is to minimize F(X 22 )/ 
the subproblem statement is as follows. 

find X 22 so that 

F(X 22 ) is minimum (3. 13. a) 


g2i(X22) ^ 0 


(3.13.b) 



(3.13.C) 



(3.13.d) 


/ • • • / 


1; V = 1, . 


• • / 


n 



(3.13.e) 


(1-6)X°2 S X 22 ^ (1+6)X22 


(3. 13.f) 


The corresponding transformed subproblem is now 


find X 


ll'^ll'^ll 


so that 


F(X 2 i ) -CE ii+drii j is minimum 


(3. 14. a) 
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gil(Xii)+£ii S 0 


-£ - - (X- T -Xt T )+Et T ^ 0 

VVll' 11 11' 11 


■'^11 ^ '^yvll^^ll'^ll^ ^ ^11 


(3. 14.b) 
(3. 14.C) 
(3 . 14.d) 


y 2/. V X/... / n 


x^i ^ Xj, s x“, 


(1-6)X°^ < X^^ < (1+6)X°^ 


^11 ^llmax 


“Ml ^ ° 


(3.14.e) 
(3. 14.f ) 
(3.14.g) 
(3 . 14. h) 


The transformation of this subproblem is almost identical to 

that of the two previously discussed subproblems. The only 

difference is that even though is still forced to 

increase, it is limited to a small positive value e^. 

This avoids overdesigning, that is, producing designs with 

excessive margin with respect to the constraints. If e^^ 

were allowed to become large and positive, the measures of 

* 

inequality constraint violation would become large 

and negative, an unnecessary requirement. The additional 
variables e^^ and are both driven to zero by addition of 

a penalty term to the objective function. The exponential 
is not required for e^^ any more as it is constrained to 
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remain negative or, at best, slightly positive. The 
positive constants c and d may be adjusted, if necessary, to 
force and within prescribed tolerances. 

3 . 5 ORGANIZATION OF THE ALGORITHM 

This section presents the organization of the algorithm. 
Within a given iteration, the general progression is top- 
down. The first iteration must be conducted without 
responses from lower level subproblems. The process 
terminates when: 1) all constraints of Problems 3.1, 3.7 and 
3.13 are satisfied, and 2) the objective function F(X^^) is 
stationary. 


i) Iteration r=l: for i from 1 to 1, for j from 1 to n^ . 

• Read initial design X° . . 

• For i>l, calculate control inputs CI^^ from control 
outputs CO a<i. 

ap 

* * * 

• Optimize subproblem ij, obtain X. ., t. n . . . 

- If i=l, solve Problem 3.14, without Eqs. 3.14.C 
and d. 

- If l<i<l, solve Problem 3.8, without Eqs. 3.8.d 
and e . 

- If i=l, solve Problem 3.3. 

• For i<l, calculate control outputs CO. . from X*.. 

ij 13 
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ii) Calculate responses: for i from 1-1 to 1, for j from 1 

to n . . 

1 

• For i>l, perform sensitivity analysis, obtain 
(Eq. 3.11 if l<i<l, Eq. 3.5 if i=l). 

• Calculate response inputs from response outputs 

RO^g, a>i (Eq. 3.12). 


iii) Iteration r=r+l: for i from 1 to 1, for j from 1 to n^ . 

• Define new initial point as previous optimum design: 

o * 

X . =x . . . 

• For i>l, calculate control inputs from control 

outputs CO a<i. 
op 


• Optimize subproblem ij, obtain X^^ 

- If i=l, solve Problem 3.14. 

- If l<i<l, solve Problem 3.8. 

- If i=l, solve Problem 3.3. 


E . . , n . . . 
ij 13 


• For i<l, calculate control outputs CO. . from X. . 

13 13 


iv) Check convergence: if converged, exit; otherwise. 


return to ii ) . 
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This organization specifies that the algorithm progresses 

top-down rather than bottom-up as suggested in Ref. 20. The 

disadvantage of conducting the optimization as proposed here 

is that the first iteration must be done without responses 

from lower level subproblems. Therefore, that iteration may 

be considered wasted as the first subproblem designs will 

not have been '-obtained with full consideration for the 

requirements from those lower level subproblems. However, 

this negative effect will be somewhat tempered if the move 

limits (Eqs. 3.8.h and 3.14.f) are kept for the subproblems 

above the lowest level. The argument for proceeding top- 

down is the need for accounting for changes in lower level 

subproblem response inputs while optimizing a subproblem on 

a given level (see term in ^*3- 3.12). If 

optimization were conducted bottom-up, the sensitivity 

derivatives making up a subproblem ROy^ would be calculated 

for a given ^^y^ would be used in a higher level 

subproblem in the same iteration before new RI values were 

yv 

generated, therefore making it impossible to calculate 

ARI . 
yv 

If the decomposition involves only two levels, the 

objection to conducting the optimization bottom-up 

disappears as there is no ARI term to calculate. But 

yv 

then, there is little difference between top-down 


and 
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bottom-up optimization. It is only a matter of deciding 
whether optimization will start at the top or at the bottom; 
from that point, the optimizations of the different levels 
follow always in the same sequence. 

3 . 6 MOVE LIMIT SELECTION 

The move limits of Problems 3.8 and 3.14 may be chosen on 
the basis of experience, or, using the bound estimates of 
Sec. 2.4. When optimization is conducted with fixed move 
limits, it is observed that the design first progresses 
steadily towards the optimum solution but does not settle; 
rather, it oscillates about that solution. During these 
oscillations, an excessive number of switches in the active 
constraint set occur, thereby preventing the process from 
converging. This situation is similar to that observed when 
solving a nonlinear mathematical programming problem as a 
sequence of linearized problems. When that type of method 
is used, similar oscillations are encountered if the optimum 
of the initial problem does not lie at a vertex of the 
feasible domain; the solutions of successive linear problems 
oscillate between adjacent vertices of that domain. As 
discussed by Pope (Ref. 27), an approach to overcome that 
difficulty is to reduce the move limits as the oscillating 
behavior appears. This approach was successfully followed 
in this study. 
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3 . 7 FINAL COMMENTS 

This chapter discusses an approach for optimizing decomposed 
design problems. It must be emphasized that no firm rule is 
established for devising a decomposition method. For a 
specific problem, the decomposition must be made on the 
basis of practical experience. As pointed out by Sobieski 
(Ref. 20) the large design problems for which this type of 
approach is conceived have always been treated in decomposed 
form. For major projects, the design process is carried out 
by the design office of a company (or even several 
companies) . That office is itself divided into specialty 
groups and subgroups each of which handles the design of a 
particular type of subsystem using its own codes and 
analytic tools. The team managing the design is then 
concerned with resolving the trade-offs. The proposed 
method offers a means for doing this with a systematic 
mathematical procedure; furthermore, it lends itself to 
parallel or distributed computer processing. 

Before applying the proposed method to a problem, the 
decomposition must be defined. This implies a definition of 
the different subproblems and the choice of their design 
variables and constraints. The structure of the controls 
must be identified as well. For the sake of generality, the 
developments of this chapter make the assumption that each 
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subproblem controls each lower level subproblem; as a 
consequence, each subproblem responds to each higher level 
subproblem. Clearly, this situation is undesirable from an 
economic standpoint. An expensive item in the scheme 
remains the calculation of the sensitivity derivatives, as 
these require calculation of second-order behavioral 
derivatives (Chap. 2). Therefore, as a general rule, the 
number of controls should be kept as low as possible. 
Practical experience should help in determining what control 
or response must be accounted for. In particular, controls 
or responses yielding typically low sensitivity derivatives 
(3E^j/3CI^j and in Eq. 3.10) should be 

ignored. 

A key point of this discussion is the approach chosen for 
finding individual subproblem designs (solution of Problems 
3.1, 3.7, and 3.13). Two additional variables are 

introduced in each subproblem. These monitor the most 
violated equality and inequality constraints. This 

particular formulation was used because i) it permits 
equality constraints to be treated numerically, and ii) it 
is compatible with the use of a usable-feasible direction 
algorithm for solving the individual subproblems. It must 
be noted that other approaches could be used as well. 
Reference 21 discusses an example based on a penalty 
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function formulation, where a cumulative measure involving 
all the constraints of the problem (except side constraints 


and move limits) is minimized. 


Chapter IV 


APPLICATION OF THE ALGORITHM TO THE DESIGN OF A 

SAILPLANE WING 

This chapter demonstrates the applicability of the proposed 
algorithm to the problem of designing a high-performance 
sailplane wing. To achieve this goal, we describe the 
proposed decomposition. We identify the subproblems, their 
design variables, and the different control and response 
vectors. The objective of the overall problem is taken to 
be the optimization of a classical performance index called 
the cross-country speed. Clearly, to efficiently design an 
airplane, one must optimize its entire configuration all at 
once. However, in order to reduce the problem size, we 
focus our efforts on the wing design, in effect solving the 
problem of redesigning the wing of an existing sailplane. 
We assume that the description of the other sailplane 
components is entirely known. The variables of the overall 
problem describe the wing shape and its structure. The 
decomposition has four levels (see Fig. 4). 

The top level of the hierarchy is occupied by a 
subproblem that selects global performance parameters. This 
drives an aerodynamic subproblem on the second level where 
the wing shape is designed to the specified performance 
characteristics. The synthesis of the wing structure is 
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Figure 4; 


Four-level hierarchy for performance sailplane 
wing design problem. 
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done on the two lower levels. The third level deals with 

the determination of spanwise distributions of global 

structural characteristics, and the fourth level is 

concerned with the detailed design of wing elements. The 

subproblems are identified by the following subscripts: p 

for performance subproblem, a for aerodynamic subproblem, g 

th 

for global structural subproblem, and Ij for g l^ocal 
structural subproblem. The developments were carried out to 
a point where computer implementation is possible. The 
details of the formulation are given in Apps . A-D. 

4.1 PERFORMANCE SUBPROBLEM 

A performance subproblem occupies the highest level of the 
hierarchy, and it drives the whole design process. Its task 
is to select a combination of performance parameters which 
maximize the sailplane’s cross-country speed. The cross- 
country speed is defined as the average speed the sailplane 
achieves over a typical flight segment. It is shown in App. 
A that the variables entering the calculation of that 
performance index are: the sailplane gross weight (W), its 

planform area (S), the coefficients (Cj^q/ *“D2 ^ ^ 

performance curve called the drag polar, and the maximum 
lift coefficient • These variables will be taken as 

the performance variables. 



{W,S,CpO'Sl''^D2'^Lmax 


(4.1) 
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The design will be subjected to the constraint that the 
sailplane be capable of a minimum rate of climb in a typical 
thermal updraft. The variables may be divided into two 
groups. On the one hand, the weight selected on the basis 
of performance considerations controls the structural design 
by determining the amount of structural material available. 
On the other hand, the remaining variables (referred to as 
aerodynamic performance variables ) control the aerodynamic 
design by specifying the aerodynamic characteristics 
required from the wing. This is, of course, a statement 
that optimization of an airplane entails carefully 
integrated aerodynamic and structural designs. The control 
outputs of the performance subproblem are therefore the 
performance variables themselves. 


CO = X 
P P 


(4.2) 


4.2 AERODYNAMIC SUBPROBLEM 


The wing 

is 

narrow, unswept, and 

twisted, and it 

has a 

double taper 

( see 

Fig. 14). Its 

airfoil is chosen 

from a 

standard 

catalog 

which details 

its parameters. 

The 


objective of the aerodynamic subproblem is the selection of 
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the wing geometry which corresponds to the values of 
planform area, drag polar coefficients and maximum lift 
coefficient specified in the performance subproblem (these 
will be referred to as target aerodynamic performance 
variables). As detailed in App. B, the aerodynamic 
variables include the wing semispan (b) and the spanwise 
position of the taper change (hj^) • Also included among the 
variables are the section chord at the root, break, and tip 
(Cr/Cb/Ct), the section maximum thickness at the root, 
break, and tip (t^,tj^,t^), and, finally, the wing twist 
angle at the break and the tip (Ej^,£^). Hence, 


X = fb,b, ,c ,c. ,c. ,t , t, , t . , e, , £ . ] 
a ‘ b' r' b t r' b t' b t* 


( 4 . 3 ) 


Once these variables are known, it is possible to calculate 
the resulting wing aerodynamic performance variables. By 
requiring that these must equal the target values specified 
in the performance subproblem, we obtain six equality 
constraints among the components of X . Furthermore, one 

3 . 

inequality constraint is added which forces the initiation 
of stall at an inboard wing section. As explained in App. 
B, the calculation of sailplane aerodynamic performance 
characteristics takes into account the effect of the 
Reynolds number, and therefore, the sailplane weight must be 
known. The vector of control inputs into the aerodynamic 
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subproblem includes then the target values for the 
aerodynamic performance characteristics and also the 
sailplane weight. 


Cl 


a 



(4.4) 


The aerodynamic design controls the structural design 
through two major elements: the wing shape and the 
aerodynamic loads. The former determines the root bending 
moment as well as the size of the load bearing part of the 
section; the latter clearly influences the total 
distribution of loads. All the aerodynamic variables are 
required for calculating the wing shape and the aerodynamic 
loads, therefore 


C°a = <^-5) 

Finally, the vector of response outputs from the aerodynamic 
subproblem is given by Eq. 3.11. 



3c 


3c 


3c 


3X 


3RI 


3X 




3RI 

a 


* \ 



(4.6) 


4.3 GLOBAL STRUCTURAL SUBPROBLEM 

The sailplane gross weight was specified in the performance 
subproblem. Since the weight of the other sailplane 
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components is fixed, the weight allocated to the wing (the 
target weight) is known. At this level, the wing is modeled 
as a straight beam. This subproblem objective is the 
selection of global wing structural characteristics, that 
is: the spanwise distributions of weight per unit length, 
bending and torsional stiffnesses. These are written as 
superpositions of interpolating polynomials, the 
coefficients of which are taken as the subproblem design 
variables and are referred to as wing stiffnesses (see App. 
C), 


X = {w ,...,w ,EI El ,GJ ,...,GJ 

g ' o 'no n o n 

s s s 


( 4 . 7 ) 


The compatibility between w^ , EI^ and GJ^ is enforced at the 
lowest level when the detailed synthesis forces the element 
designs to match the global level stiffnesses. The design 
is subjected to one equality constraint that specifies that 
the total wing weight equals the target weight selected in 
the performance subproblem. In addition, two global 
structural inequality constraints are added. The wing 
bending displacements must be limited, and its torsional 
divergence speed must exceed a given minimum. The 


calculation of these constraints requires the wing shape, a 
number of aerodynamic load parameters and, of course, the 
target weight. The wing shape (vector H) is a function of 
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the aerodynamic variables, while the target weight (W) is a 
function of the performance variables. The aerodynamic load 
parameters (vector A) are themselves functions of the wing 
shape and the sailplane weight as shown in App. B. Hence, 
the vector of control inputs into the subproblem is 


CI^ = {W(X ),h'^(X ),a'^(X ,W(X ))} 
g ' ' p' ' ' a^ ' a ' p' ' ’ 


(4.8) 


The choice of the wing stiffnesses controls the design of 
the details of structure, hence 


CO = X 

g g 


(4.9) 


The response outputs from the global structural subproblem 
are given by Eq. 3.11. 



f 

k 

3e 

k 

3£ 

* 

as 

* 

3e 

* 

d£ 
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g 

g 
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^ ' 

3W 

3H 

3A 

3RI 

3X° 
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* 
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* 

g 

g 

g 

g 

g 


(4.10) 


^ ' 

3W 

3H 

3A 

3RI 

3X° 






g 

g J 



4 . 4 LOCAL 

STRUCTURAL : 

SUBPROBLEMS 




The fourth 

level 

of the hierarchy 

is devoted 

to designing 

the details 

of the 

wing 

to 1 

match the 

stiffness 

distributions 


chosen by the global level subproblem. To achieve this, the 


wing is divided into n constant property elements (see Fig. 

s 
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22), each of which is designed in a separate subproblem. 
The design variables describing one element include the 
cross-sectional dimensions of the spar caps (a,t^), the 
thicknesses of the sandwich panels constituting the leading 
and trailing edge shells spar webs 

(t^,t^^), and the rib spacing (d). Hence, 


XT. = |a.,t .,t . , t,- .,t .,t^ .,d. 

I3 J cj ' sj ' fsj' wj fwj j 


(4.11) 


The variables are subjected to three equality constraints 
insuring that the element stiffnesses equal the target 
values (vector ) specified at the global structural level. 
Furthermore, it is required that the elements withstand 
applied end loads (vector ) without failure or local 
buckling. App. D shows that the calculation of the 
stiffnesses involves the element gross dimensions (vector 
) ; these are functions of the aerodynamic variables. 
Also, the applied end loads (M^) computation must be based 
on the vector describing the wing shape (H), a function of 
the aerodynamic variables, the vector of aerodynamic load 
parameters (A), a function of aerodynamic variables and 
weight, and the global structural variables (Xg) to account 
for the effect of inertia loads. Hence, 


Cl 


1] 


{k'^(X ),h'^(X ),m'^[X ,A(X ,W(X )),H(X )]} 
3 3 ^ 3 g a' ' p' ' ' a' ’ 


(4.12) 
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Equation 3.5 gives the response outputs from subproblem 1 j , 
a lowest level subproblem. 


RO 


Ij 


3e 


Ij 


•Ij' 


3K . 
J 


* * 

d £ ^ d £ ^ 

l3 Ij 


, T) 


3H . 
3 


3M 


13 ' 


3,*. 


3K . 
3 


3tii. 3n^. 


3H . 
3 


3M 


(4.13) 


4.5 RESPONSE INPUTS TO THE DIFFERENT SUBPROBLEMS 
The vector of response inputs into a subproblem above the 
lowest level is given by Eq. 3.12. For the global 
structural subproblem, we have 


T * T* * T * T 

^^g ^ ^^^llg^ ' • • • ' ^^^In^g^ ' ^’’llg^ '• 


s^ 


where, for example 


* T 

[t, . } = 

i3g 


< 


* 

T, z, . 
Ig Ij 


3e* 3M. 3e* 3K. 

l3 3 l3 3 


3M.3X 

3 g 


3K .3X 

3 g 


(4.14) 


The vector of response inputs into the aerodynamic 
subproblem is 


RI^ = { {e 1^, {e, , j^, . . . , | e, 1^, 

a ' ' ga‘ ' ' lla‘ ' ' ' In a’ ' 

s 


* T * T * T 

^^ga^ ' ^’’lla^ / • • • / 1 

s 
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* T 

1 = 
‘ ga ‘ 


/. 3e 3e 

fs + — ^ARI + — ^AX > , 
gal g g ^ g 




3RI 


3X 


3e 3A 3e 3S 

__3 / _g 

3A 3X 3S 3X 

a a 


and 


I'ljal = 


* 

'I'l El • 
la Ij 


3e, -3M .3A 
Ij J 


3M . 3A3X 
J ' a 


3e, .3M.3H 3e, .3H. 

Ij J Ij J 


3M.3H 3X 
J a 


3H.3X 
J a 


(4.15) 


Finally, the response inputs into the performance subproblem 
are 


RI^ = {{e {e r,{E,T }™, (e, 

p ‘‘ ap‘ gp ‘lip’ * 

*T*T*T * T 

|T1,„ 


ap- 


gp 


In p' 
s^ 


* T 

{e ] = ^e + 

* ap’ ) a 


* * 

. 3e 3e 

^ ARI + ^AX° 

a a 


3RI 


3X 


3e*3W 

a 


3W 3X 


Also 
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* T 

gp 


^'1'^ (e* + 

gp I g 


u 


ARI + 

g 


8RI 



as* 3A 

g 

3A 3X 


3e 3W 

g 


3W 3X 


and, finally 


|e, . } ■" = ly, t, . 

1JP‘ ^ Ip IJ / 


* 

3M . 3A3Wl 
3 \ 

3M . 

3A3W3X 1 

3 

Pj 

.4.d 

and 3 . 1 


(4.16) 


According to Eqs . 3.4.d and S.lO.d., we must further require 


4', +t, = 1 

Ig la Ip 


T +T = 1 
ga gp 


(4.17) 


Note that we have specified = 1 (Eq. 4.16) as there is 
only one subproblem controlling the aerodynamic subproblem. 
Furthermore, we assume the participation of each subproblem 
into the local structural subproblem designs to be equal, 
for example 


4' =4' 

llg 12g 


= 4-, = 4't 

In^g Ig 


The coefficients 4"^ and 4', are defined similarly 

la Ip ^ 


66 


s 


CRiOJf^AL' PASS ro 
OF POOR QU.ALrn' 


3 . 1 FINAL REMARK 

This decomposition conforms to the general form defined in 
Chap. Ill; in particular, it does not contain reverse or 
network controls. This is because the following effects are 
neglected (see Apps. B and C) : i) effect of wing 
flexibility on aerodynamic characteristics, ii) effect of 
changes in wing section elastic axis position on torsional 
divergence speed, and, iii) effect of changes in wing 
section center of gravity position on inertia-force-induced 
torque. To account for i), one must include reverse 
controls between the global structural and aerodynamic 
subproblems. Effects ii) and iii) require reverse controls 
between local structural and global structural subproblems. 


0R1G'?^AI.' IS 
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Chapter V 

TWO-LEVEL MINIMUM WEIGHT WING DESIGN 

A complete numerical application of the algorithm presented 

in Chap. Ill is appropriate, but it was decided that 

computer implementation of the four- level high-performance 

sailplane design problem formulated in Chap. IV and Apps. 

A-D would be a task well beyond the scope of this research. 

A more limited and manageable but still significant 

application of the algorithm is the minimum weight design of 

the sailplane wing. This corresponds to the two lower 

levels of the hierarchy described by Fig. 4, except that the 

objective function of the overall problem, and therefore of 

the top level of the decomposed problem, is now the total 

wing weight. The global structural subproblem now has the 

form of Problem 3.13. The objective function (F(X^^) = 

)) is the weight, and there are three inequality 

constraints (gii(Xii ) = gg(Xg)); two set upper and lower 

limits on the wing tip displacement, one provides a lower 

bound for the wing torsional divergence speed (see App. C). 

The weight and stiffness distributions are Lagrangian 

interpolations defined on five equal length intervals (n = 

s 

5), thereby requiring six coefficients per distribution, for 
a total of eighteen global level variables (X ) . The five 

g 
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local structural subproblems are formulated as in Problem 
3.1. Each one is described by seven design variables ) . 
The vector of constraints contains twelve 
constraints on stresses, five on buckling (only the buckling 
of the sandwich panels is considered, that of the spar caps 
is ignored, see App. D) . Also three constraints are 
included which limit the ratios of spar cap width to element 
chord and of sandwich panel total thickness to face sheet 
thickness (see App. E). Finally, three equality constraints 
(ej^^(X^j)) require that the element weight per unit length, 
bending and torsional stiffnesses equal the target values 
specified at the global structural level (see App. D) . 

The multilevel decomposition for this two-level problem 
corresponds to the two lower levels of the hierarchy of Fig. 
4 except that the definition of the wing planform is now 
assumed fixed so that the vectors of wing shape (H) and 
element gross shape (H^ ) are constant throughout the 
optimization. Furthermore, the effect of the sailplane 
weight (W(Xp)) on the aerodynamic parameters is ignored so 
that vector A is constant as well. The control inputs into 
local structural subproblem Ij (CI^^) are then the vectors 
of target stiffnesses (K^ ) and end loads (M^ ) . The response 
outputs from that subproblem are given by Eq. 4.13, without 
the sensitivity derivative with respect to . The response 
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inputs into the global structural subproblem are given by 
Eq. 4.14 where the participation factor is set to unity 
as the global structural subproblem is the only one that 
controls the local structural level subproblems. The 
numerical values of the parameters needed for the 
calculations of Apps. A-D are given in App. E. 

5 . 1 TEST DESCRIPTION 

Two test cases will be described. The first starts from the 
feasible domain, all the constraints being satisfied with 
ample margins. The second starts from the infeasible 
domain. It is such that, at the global level, the tip 
displacement exceeds its upper limit by 130%, while the 
divergence dynamic pressure is 10% below its lower limit. 
With the exception of the tip element, all the wing elements 
have a large proportion of stress and buckling constraints 
violated, resulting in measures of inequality constraint 
violation ranging from -3.0 to -6.0. In both cases, 
an initial design was given for each wing element and the 
stiffnesses of these elements were calculated. Trial and 
error was then used to determine the global level 
stiffnesses whose averages would yield back the local level 
stiffnesses just calculated, in accordance with Eq. C.23. 
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Figure 5: 


Two- level wing design, weight convergence 
history. 


71 


Figure 5 gives the convergence history for the test 
cases. Both optimizations took 28 iterations^ to produce a 
feasible design with a stabilized weight of 1030 N. (for two 
wings), with the differences in the final values being less 
than 1%. The design started in the feasible domain achieves 
a weight reduction of 68%. That begun in the infeasible 
region gains barely 6% of weight as it redistributes the 
structural material in an attempt to satisfy the violated 
constraints. Both final solutions have the divergence 
constraint critical, while the upper displacement constraint 
is satisfied with a margin of about 25%. The designs of 
wing elements 1 (inboard), 2 and 4 are critical or near 
critical (-.01 < e S .07), while those of elements 3 and 5 
(outboard) are not (t > .10). The active constraints for 
the critical elements are: i) the stress constraints in the 
carbon fiber spar caps, and primarily in the lower one, ii) 
the shear buckling constraint for the upper and lower 
sandwich panels of the trailing edge shell and the rear spar 
web, iii) the constraint setting an upper limit on the spar 
cap width to section chord ratio (a/c). Furthermore, for 
each element, the rib spacing d is against its upper limit. 


^ Unless otherwise specified, one iteration corresponds to 
one complete optimization of the overall problem. It 
entails one optimization of the global structural 
subproblem, one optimization and one sensitivity analysis 
for each local structural subproblem. 
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while the spar web face sheet thickness is against its lower 
limit. The parameter b in the local level subproblem 
objective function (Problem 3.3) was set to 100. At the 
global level, the objective function parameters (c and d, 
see Problem 3.14) were chosen so that, at the beginning of 
each global level optimization, the penalty term in the 
objective function ( -ce^^+dii^^ ) was ten times the weight 
(F(X^^)). As a consequence, all equality constraints were 
satisfied, after most system iteration, to within 1%, or 
less { n S . 01 ) . 

5.2 COMPARISON OF GLOBAL AND LOCAL DESIGNS 

Figure 6 shows the global and target local weights for the 
initial and final designs of both tests. The same 
information is given for the bending stiffnesses on Fig. 7. 
On these figures, the straight-line segments join 
stiffnesses corresponding to the same (initial/final) 
design; they do not represent the global level stiffness 
distributions, as these are obtained with Lagrangian 
interpolations and have no slope discontinuities. The 
figures show good agreement between the final target values 
(open symbols) toward the wing root, and the agreement 
worsens somewhat toward the wing tip. The variability is 
even more marked with the final global level variables 


73 


Weight per 
unit length 
(10^ N/m) 


Figure 6: 


ORlGif^IAL PAGE rs 
OF POOR Q’JAijjY 


2 . 5 ^ 



Span (m) 


Legend : 


• 

O 

A 

A 


Initial values 
Final values. 
Initial point 
Initial point 
Initial point 
Initial point 


feasible; global weight, 
feasible; element target weight, 
infeasible; global weight, 
infeasible; element target weight. 


Two-level wing design, comparison of initial and 
final global and target weights per unit length. 
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Bending 
stiffness 
(10^ Nm2) 
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Legend: Initial values. 

Final values. 

• Initial point feasible; global bending stiffness. 

O Initial point feasible; element target bending stiffness. 

▲ Initial point infeasible; global bending stiffness. 

A Initial point infeasible; element target bending stiffness. 


Figure 7: 


Two-level wing design, 
final global and target 


comparison of initial 
bending stiffnesses. 


and 
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(solid symbols). This is probably evidence that the optimum 
is fairly "flat". In other words, the global level design 
may be perturbed without drastically changing the overall 
design quality measured in terms of minimum objective 
function and constraint satisfaction. Three factors may be 
contributing to that situation. First, due to the 
relationship between the number of coefficients in the 
global level wing model (6 per distribution of weight per 
unit length, bending and torsional stiffnesses) and the 
number of wing elements (5), there is an infinite number of 
global level designs which yield the same local level target 
values. Second, it may be proven algebraically that all the 
combinations of global level weights which result in the 
same local level target weights actually result in the same 
total wing weight (this can be shown easily for Lagrangian 
interpolations based on an odd number of equal length 
intervals). Third, all the global level constraints (in the 
untransformed subproblem) are integral constraints (see App. 
C) . The integration process should be expected to smooth 
out the effect of variations in global level stiffness 
distributions . 

The distributions of stiffnesses are wavy; this must be 
related to the order of the interpolating polynomials used 
to model them. These polynomials are of the fifth order in 
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the present model. It would be interesting to model the 
wing weight and stiffness distributions with lower-order 
polynomials; this should reduce the wavyness while 
necessitating fewer variables in the global structural 
subproblem and making the global and local structural models 
independent. This would, however, require a new definition 
of the relationship between the global level weights and 
stiffnesses and the corresponding local level target values 
(Eq. C.23). 

Table 1 compares numerically the initial and final 
designs for the inboard and outboard elements. As expected, 
the trend observed among the target stiffnesses is repeated 
at the lower level; the two designs obtained for the inboard 
wing element are very close, while those for the outboard 
element differ substantially. This general trend can be 
intuitively understood as the overall wing response must be 
more sensitive to the design of the wing root which carries 
the total wing load than to that of the wing tip which is 
almost unloaded. 

5.3 QUALITY OF SENSITIVITY-DERIVATIVE-BASED EXTRAPOLATIONS 
Figure 8 depicts a typical evolution of optimum measures of 
inequality and equality constraint violation and ^ • 
The example corresponds to the design of the second wing 


TABLE 


Two-level wing design, comparison of initial and final 
designs for the inboard and outboard elements. 


Wing 

Element 

I nboard 

Outboard 

Initial Point 

Feasible 

Infeasible 

Feasible 

Infeasible 

Design 

Initial 

Final 

Initial 

Final 

Initial 

Final 

Initial 

Final 

★ 

a 

xlO'^ 

3.00 

9.01 

3.00 

9.22 

3.00 

3 . 44 

3.00 

1.28 

t 

c 

xlO- 3 

20.00 

2.48 

3.00 

2.37 

20.00 

4.38 

3.00 

11. 10 

^s 

xlO" 3 

10.00 

5.98 

3.00 

5.96 

10.00 

3.92 

3.00 

3.91 

^fs 

xlO"* 

20.00 

5. 13 

5.00 

5. 15 

20.00 

7.13 

5.00 

5.22 


xlO'3 

7.00 

9. 57 

3.50 

9.50 

7.00 

2.76 

3 . 50 

7.61 

^fw 

xlO"* 

10.00 

* * 

5.00 

5.00 

* * 

5.00 

10.00 

* * 

5.00 

5.00 

* * 

5.00 

d 

xl0‘ ' 

2 . 50 

* * 

5.00 

5.00 

★ * 

5.00 

2 . 50 

* ★ 

5.00 

5.00 

*• -k 

5.00 


All dimensions are in meters. 


Variables against side constraints. 


OF POOR QUALiTY 


78 


ORiSlNAL PAGS IS 

OF. POOR QUALITY 



I terations 


Two-level wing design, typical history of 
measures of constraint violation. 

Wing element 2 , infeasible initial design. 


Figure 8: 
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element when optimization is started from the infeasible 

domain. A short digression is necessary to explain the 

physical meaning of variables e and r\. For an infeasible 

design, e measures the violation of the most violated 

constraint; for a feasible design it gives an indication of 

the margin associated with the constraint which is the 

closest to being critical. For example, if we refer to a 

design problem involving a single inequality constraint on 

displacement u (g(X) = u(X)/u -1^0), a value of -1.5 for 

E indicates that displacement u exceeds its maximum value 

u by 150%, while a value of 0.5 means that the 

m 3.^ 

displacement is only 50% of its upper limit. Variable t) 
measures the violation of the most violated equality 
constraint of a subproblem. Considering a subproblem 
subjected to the requirement that a function f of its design 
variables remain equal to a target value f^ (e(X) = 
f(X)/f^-l =0), a value of .25 for t) means that f(X) is 25% 
away from f^. 

Returning to Fig. 8, the solid lines give the evolution 

•k k 


of e ^2 ^12' T'hey describe the evolution of these 

optimum quantities throughout the design process, each value 
being obtained at the end of a complete optimization of 
subproblem 12. The dots correspond to the estimates of 
these as generated on the basis of optimum sensitivity 
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analysis in the global structural subproblem. They are, 

+ * ★ ★ 

respectively, e = e^- and ti .. = in the 

^ yvll 12g vvll 12g 

notation of Chap. Ill (see Problem 3.13). The 

interpretation of the graph may be aided by an example. At 

* 

the end of iteration 3, e^^ is -1.35; the global 

structural subproblem improves its design at the beginning 

of iteration 4 and predicts that this will result in a 

* 

change to -.68 while e^^ actually is reduced to -.78 when 

local structural subproblem 12 is optimized. 

From Fig. 8, one can observe that the algorithm 

effectively reduces the violation of inequality constraints 

while keeping the violation of equality constraints within 

1%, most of the time. The success is a result of the 

generally good prediction of these measures by the global 

structural subproblem. Occasionally, however, the 

predictions are poor, as for iterations 6, 22 and 24 for 

* * 

e^ 2 ' '’^12' predictions can be 

correlated with changes in the set of constraints 
determining the successive solutions to subproblem 12; as 
discussed in Sec. 2.4, these changes do introduce 
discontinuities in the sensitivity derivatives. This slows 
down the convergence process, but the algorithm appears 
robust enough to overcome the difficulty; this observation 


QUftUTY 


is consistent with the results of Ref. 21. It must be noted 
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that the occurrence of active constraint switching does not 
always produce discontinuities as significant as those 
observed here. As a matter of fact, about 60% of the 
optimizations included in the example presented end with 
changes in the active constraint set, but only 15% seem 
really affected by the change. 


5.4 ADDITIONAL COMMENTS 

As was pointed out in Chap. Ill, limits must be set on the 
global level variables to insure that the sensitivity- 
derivative-based linear extrapolations of the local level 
measures of constraint violation are not used beyond their 
range of validity. Upper and lower bounds could be 
estimated for these limits, using Eqs. 2.25 and 2.16. 
However, it was elected not to do so in this study. The two 
test cases presented here were started with move limits 
allowing changes of ± 10% in the global level variables at 
each iteration of the overall problem (6 = .10 in Eq. 
3.13.f). As the process neared the optimum solution, these 
limits were reduced as discussed in Sec. 3.6. For example, 
in the test case started from the infeasible starting point, 
the move limits were kept at 10% during the first 5 
iterations, 5% for iterations 6-21, 2.5% for iterations 
22-25, and finally, 1% for the last three iterations. It 
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must be pointed out that no attempt was made to refine the 
bound selection strategy nor to identify the widest usable 
bounds for the first steps of the process. Doing so could 
probably improve the convergence characteristics of the 
algorithm. 

The impact of active constraint switching on the rate of 
convergence is related to the specifics of the solution of 
the individual subproblems. This approach tracks the most 
violated equality and inequality constraints. Therefore, a 
relatively large number of changes in the active constraint 
set may be expected. By contrast, the approach of Ref. 21 
combines most constraints (except upper and lower bounds and 
move limits) in a penalty function which is continuous in 
the constraints. That approach should not result in so many 
changes. The question of which approach provides the best 
rate of convergence may be answered only by using both 
approaches to solve the same design problem. 

Early tests with this algorithm revealed a situation that 
deserves comment. In some circumstances, the subproblem 
solution may not be unique. An example of this is depicted 
on Fig. 9 where the subproblem involves one equality 
constraint that is parallel to one of the sides of the 
domain defined by the inequality constraints. Assuming that 
the equality constraint is satisfied at the outset of the 
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Figure 9: 


Situation where a subproblem has an infinite 
number of solutions. 
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optimization, we see that each point of segment AB is a 
solution. This situation does not cause difficulties 
during optimization of the subproblem. However, it leads to 
singular sensitivity matrices. This can be shown on the 
example. As none of the subproblem objective function 
(e ^+bTi) or active constraints depend on variables the 

corresponding row and column of the sensitivity equations 
(Eq. 2.14) are identically zero. This difficulty may be 
overcome by eliminating the null row and column from the 
equation. In this problem, the situation appeared when the 
design process was started with all variables at lower 
bounds . 

To conclude this chapter, a few words about the 
implementation are in order. A special purpose FORTRAN 
computer program was developed to carry out the examples 
described. It was run in batch on an IBM 3081 processor. 
Three iterations were run at a time and results from 
successive iterations were kept on a direct access data file 
so as to provide the capability for restarting the 
optimization from any iteration. Individual subproblems 
were optimized with the usable-feasible direction algorithm 
implemented in subroutine CONMIN (Ref. 28). For the two 
examples discussed, convergence was achieved in 28 


iterations; a typical iteration took 30-35 CPU seconds. Of 
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that time, about 60% was devoted to optimizing the global 
structural subproblem, while the remaining 40% was split 
evenly between optimization and sensitivity analysis of the 
5 local level subproblems. During one iteration of the 
overall problem, optimization of the global level subproblem 
took an average of 18 iterations to converge. Each global 
level iteration required about 21 function evaluations 
(weight, wing tip displacement and wing divergence dynamic 
pressure). At the local level, a wing element design took 
an average of 38 iterations. Each local level iteration 
required about 11 function evaluations (stiffnesses, 
stresses and buckling constraints), while each local level 
sensitivity analysis required 531 function evaluations. The 
numbers of function evaluations reported above include those 
required for derivative calculations. For reasons of 
convenience, all derivatives for optimization were obtained 
by forward difference formulae; those for sensitivity 
analysis were found by central difference formulae. 
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Chapter VI 

DISCUSSION AND RECOMMENDATIONS 

This work centers on a multilevel decomposition approach to 
the optimization of large problems. The starting point is 
Sobieski's original idea (Ref. 20) of using optimum 

sensitivity analysis to account for the coupling between the 
different subproblems of a decomposed problem. The 
objective of this research was the investigation of the 
applicability of Sobieski's algorithm to the design of 
complex engineering systems. The three major contributions 
of this dissertation are: i) a presentation of an extended 

version of the algorithm, ii) an illustration of the 
applicability of the algorithm to a realistic 

multidisciplinary design problem, and, iii) an application 
of the proposed algorithm. This chapter summarizes these 
developments and offers some recommendations. 

6.1 ON THE PROPOSED ALGORITHM 

The algorithm we developed is very general; in contrast with 
previously introduced multilevel approaches, it does not 
make restrictive assumptions as to the mathematical nature 
of the problem considered nor is it developed for a specific 
class of physical problems. Rather, it is devised for the 
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optimization of large systems which may be decomposed into 
subsystems . 

We hypothesize an existing decomposition of the overall 
problem into smaller subproblems. It is assumed that the 
decomposition is multilevel and hierarchical; that is, the 
subproblems are organized on levels corresponding to 
progressively more refined models of the system designed 
and, furthermore, the subproblems on a given level exert 
control on the design of the subproblems on lower levels. 
The controls must be defined with the hierarchy. It is 
finally assumed that the subproblems entail finding a design 
that satisfies a set of equality and inequality constraints. 

An iterative procedure is then used where each subproblem 
is designed separately by minimization of a measure of 
constraint violation, for fixed values of the controls 
imposed by higher level subproblems. Sensitivity 
derivatives of the resulting subproblem solution are 
generated which constitute the response of that subproblem 
to the higher level subproblem controls. Each subproblem is 
designed with additional constraints which require 
improvement of the lower level subproblem responses. This 
algorithm specifically allows for the extreme situation 
where a subproblem controls each lower level subproblem, 
and, as a consequence, each subproblem responds to each 
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higher level subproblem. The highest level of the hierarchy 
is occupied by a single subproblem which drives the whole 
design process. Its task is the optimization of the overall 
problem objective with respect to variables which describe 
the system in very general terms. This subproblem also 
manipulates its controls to improve the lower level 
subproblem designs. 

The specifics of the solution of the individual 
subproblems differ from those of Ref. 21, where equality 
constraints are handled analytically by elimination of 
dependent variables, while the independent variables are 
found by minimization of a penalty function. The present 
approach introduces two additional variables which measure 
i) the violation of the most violated inequality constraint 
(or the margin of the constraints which are the closest to 
being critical), and ii) the violation of the most violated 
equality constraint. These variables are then combined in a 
cumulative measure of violation which is minimized to find 
the "best" subproblem design. The advantage of this 
approach is that it is more general in that it handles 
equality constraints numerically. Furthermore, it is 
compatible with the use of a usable-feasible direction 
algorithm. 
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This discussion shows that when the decomposition 
contains more than two levels, it becomes necessary to 
consider the sensitivity of a given subproblem above the 
lowest level not only to the control inputs into that 
subproblem, as shown in Ref. 20, but also to the response 

inputs into that subproblem, as well as to the initial 
design for that subproblem. Of course, it may be argued 
that ignoring such effects will only slow down the 
convergence of the overall process while resulting in fewer 
sensitivity analyses and in an easier implementation. This 
question must be resolved by testing a simple three-level 
problem for which a solution is known or easy to generate. 

All that is required before using this algorithm is a 
choice of the decomposition of the overall problem into 

subproblems, and of the controls between the subproblems. 
These choices must be made on the basis of practical 

experience accumulated with the type of systems considered. 
The decomposition must correspond to the general description 
given above. Intuitively, the number of controls should be 
kept to a minimum; the decision of whether to include a 

specific control or not could be helped by analysis of the 
sensitivity of the optimum of the subproblem being 
controlled to the control considered. 
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6.2 ON THE APPLICABILITY OF THE ALGORITHM 

The applicability of the algorithm was illustrated by 
formulation of a realistic multidisciplinary design problem: 
the simultaneous aerodynamic and structural design of 
sailplane wing for maximum cross-country speed. A 
multilevel decomposition conforming to the desired general 
form was developed. The highest level is occupied by a 
subproblem which selects the sailplane performance 
parameters that maximize the cross-country speed. Then 
follows an aerodynamic subproblem which determines the wing 
shape that produces the aerodynamic performance parameters 
selected at the level above and meets a stall-related 
constraint. On the third level, a global structural 
subproblem chooses the distributions of wing weight and 
stiffnesses which satisfy constraints on wing tip 

displacement and divergence dynamic pressure for a total 
weight equal to that selected in the performance subproblem. 
Finally, at the lowest level, a number of structural 
subproblems are concerned with finding wing element designs 
which yield the stiffnesses specified by the global 
structural subproblem while satisfying stress and buckling 
constraints . 

The controls of this decomposition were chosen by 
considering the rigorous functional dependence of a given 
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subproblem optimum on the variables of the higher level 
subproblems. As a consequence, the weight, a performance 
variable, controls all subproblems below the highest level; 
likewise, the wing shape vector, a combination of 
aerodynamic variables, controls the two structural levels. 
This underscores the need for an algorithm that allows each 
subproblem to control each lower level subproblem. 

The formulation of this design problem was carried to a 
point where computer implementation is possible. Its 
completion can prove to be a very instructive exercise, 
especially as 'che initial theoretical developments related 
to the proposed algorithm show that handling problems with 
more than two levels may require calculations that were not 
necessary for two- level problems. 


6.3 ON THE EXAMPLE PROBLEM 

The proposed algorithm has been applied to the two- level 
minimum weight design of the sailplane wing. The tests show 
that it is effective at producing low weight feasible 
designs whether optimization is started from the feasible 
domain or the infeasible one. 

The numerical experiments conducted confirm that as 
optimization progresses, changes occur in the sets of 
constraints defining the optima of the different 
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subproblems. These changes occasionally cause 
discontinuities in the sensitivity derivatives, and these 
may even result in oscillations about the optimum solution. 
Although the algorithm successfully overcomes these 
discontinuities, they still result in delayed convergence. 
The susceptibility of a particular algorithm to active 
constraint switching depends on the specific approach used 
to optimize the individual subproblems. This point should 
be investigated carefully prior to undertaking any large- 
scale implementation of the algorithm. A good basis for 
comparison could be obtained by solving a small two-level 
problem using all candidate subproblem formulations. 
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Appendix A 

PERFORMANCE SUBPROBLEM 

The design of an airplane is always driven by the desire to 
obtain the best possible value of a given performance index. 
In this application, we will maximize the sailplane best 
cross-country speed, that is, maximize the best average 
speed the sailplane can achieve over a segment of flight 
which starts with a high-speed glide and follows with a 
circling climb in a thermal to recover the height lost 
during the glide. In addition, we will require that the 
sailplane achieves a minimum rate of climb in a typically 
weak thermal. The design variables optimized in this 
subproblem are the sailplane weight, wing area, maximum lift 
coefficient, and the parameters relating its drag 
coefficient to its lift coefficient. 

This appendix describes the calculations required to 
obtain the sailplane best cross-country speed and rate of 
climb as required for the optimization of the performance 
subproblem. 


A. 1 NOMENCLATURE 

Cj^ : sailplane drag coefficient. 

Cdq/ ^D2* f of the drag polar (Eq. A.l). 

: sailplane lift coefficient. 

Li 
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D 

d 


h 



S 

t. 


V 

V 

c 

V 


cc 


"^ith 


V 


V 


sc 


I 

V 

sc 



: drag force. 

: horizontal distance flown in typical flight segment 
(Subsec. A. 3. 2) 

: height lost and recovered in typical flight 

segment . 

; acceleration of gravity. 

: Eq . A . 8 . 

: lift force. 

: radius of thermal or circling flight. 

: radius where thermal strength is reduced to one 

half of the maximum strength v., 

^ thmax 

: wing planform area. 

: time expended in gliding flight. 

; time expended in climbing flight. 

: horizontal flight speed. 

: sailplane rate of climb in thermal flight. 

: cross-country speed. 

: rate of sink of air between thermals. 

: sailplane rate of sink with respect to ambient air 

in straight flight. 

: sailplane rate of sink with respect to ambient air 

in circling flight. 

: projection of v in sailplane plane of symmetry. 

s c 

: absolute flight speed. 

: rate of climb of air in thermal. 
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^thmax 

W 

<t> 

p 


maximum value of v. , . 

tn 

sailplane gross weight, 
sailplane bank angle, 
mass density of air. 


A. 2 SAILPLANE PERFORMANCE CURVES 
A. 2.1 Drag Polar 

The performance of a sailplane in gliding flight may be 
described by a curve that relates total drag coefficient to 
total lift coefficient. That curve is named the drag polar, 
it is essentially parabolic in shape; we will assume the 
following form 


= 


So‘"'^D1^l‘"^D2^L 


(A.l) 


where the coefficients C^q, and 0^2 depend mainly on the 

sailplane aerodynamic design and on its weight. 


A. 2. 2 Speed Polar 

As a sailplane flies unaccelerated on a straight trajectory 
in undisturbed air the lift and the drag combined 
equilibrate exactly the weight (see Fig. 10). The 
similarity between the triangles of speeds and forces yields 

Vg/V = D/L = (A. 2) 


or, with Eq. A.l 
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Figure 10: 


Sailplane in unaccelerated flight, definition of 
forces and speeds. 
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The equilibrium of the forces acting on the sailplane yields 
W = (l 2 .d 2 )V 2 

= . 5 pSV^(C^+C^)^/^ (A. 4 ) 

while, from the speed triangle, we know 
2 2 2 

v:^ = v"^+v; (A. 5) 

or, combining Eqs. A. 2 , A . 4 and A . 5 
W = . 5 pSV^Cj^[l+(C^/Cj^)^]^/^ (A. 6 ) 

For most flight attitudes 


>> 


and 


W = L = . 5 pSV C. 


(A. 7 ) 


The lift coefficient is thus related to the horizontal 
flight speed by 


Cl = 


K = 2 W/pS 


(A. 8 ) 
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So that the sailplane rate of sink (Eq. A. 3) is now related 
to its horizontal flight speed as follows 


This equation is the speed polar. It is quite important to 
the sailplane pilot since it allows him to know the altitude 
loss per unit time corresponding to his horizontal flight 
speed. 

A. 2. 3 Circling Polar 

A sailplane circling at constant speed in undisturbed air is 
depicted in Fig. 11. It flies on a helical trajectory of 
constant radius r about a vertical axis. In the plane of 
symmetry of the sailplane, the similarity between the 
triangles of speeds and forces yields 


In the vertical plane containing the center of gravity of 

the sailplane, we see that the vertical rate of sink v is 

sc 

now related to the rate of sink in the plane of symmetry 




(A. 9) 


v^^/V = D/L = C^/C^ 


(A. 10) 


V 


sc 


and the bank angle <t> by 



(A. 11) 


so that 
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Figure 11: 


Sailplane in constant speed circling flight, 
definition of forces and speeds. 
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(A. 12) 


Following a development similar to the one conducted in the 
previous subsection, we relate the sailplane lift to its 
weight by 

L = .SpSV^Cj^ = W/COS0 (A. 13) 


and the horizontal flight speed is now given by 
V= [K/(Cj^coS(j) (A. 14) 


where K is defined in Eq. A. 8. The vertical rate of sink is 
obtained by combining Eqs . A. 12 and A. 14 


V 

sc 




1/2 


(Cj^cos<> )^'^^ 


(A. 15) 


The bank angle 0 is found from the inertia forces 

W 

COS0 = 

[W^+(WV^/gr)^ ] 

= [l+(V^/gr)^]“^/^ (A. 16) 


We may eliminate V from Eqs. A. 14 and A. 16 and that yields 
COS0 = [ 1- (K/Cj^gr)^ ] 


(A. 17) 
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Finally, from Eqs. A.l, A. 15 and A. 17 we obtain 

V3^ = (A. 18) 

[C^-(K/gr)^]^/^ 

Equation A. 18 is known as the circling polar. For a given 
sailplane, it relates the sailplane rate of sink in circling 
flight to the radius of turn and the specific lift 
coefficient selected by the pilot. 

A. 3 A SAILPLANE PERFORMANCE INDEX; THE CROSS-COUNTRY SPEED 
A. 3 . 1 Soaring 

A sailplane is usually towed to altitude by another 
airplane, a winch or even a car. If the pilot encounters 
only motionless air, he will trade his altitude for the 

speed needed to sustain flight. To keep his altitude and 

really soar, the pilot must find regions of rising air. 
Three soaring techniques are frequently used: thermal 

soaring, ridge soaring and wave soaring (See Ref. 29). 

In thermal soaring, the source of lift (lift is meant 
here as upward force) is columns of rising air. During 

sunny days, radiation from the sun heats the surface of the 
earth which, in turn, heats the air in contact with it. The 
heating of the air is conditioned by the nature of the 

underlying surface. By buoyancy, the warmer air rises with 
respect to the cooler air. The regions of rising air are 
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called thermals. Once a pilot finds a thermal, he usually 
tries to circle in it to gain altitude. 

When the wind blows on a mountain ridge, it is deflected 
upwards. By flying in a direction generally parallel to the 
ridge, on its upwind side, the pilot can remain at about 
constant altitude for as long as the ridge is uninterrupted, 
provided it is reasonably straight. This is called ridge 
soaring. The maximum altitude achievable in ridge lift is 
only a few hundred meters over the mountain. 

When a mass of air flows across a mountain top, a system 
of standing waves may develop downwind from the ridge. 
Regions of rising and descending air alternate. Those 
regions may extend very high in altitude, up to ten times 
the height of the ridge. This source of lift is always used 
to establish altitude records. This type of flying is 
referred to as wave soaring. 

From this very brief discussion of soaring, it is clear 
that a sailplane may be flown at quite different flight 
regimes. Therefore, it is difficult to define a performance 
index for the design of a sailplane without making some 
assumptions as to the type of soaring it will be used for. 
In this work, we will design the sailplane for thermal 


soaring . 
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A. 3. 2 The Cross-Country Speed 

To cover a given distance, a pilot alternates slow circling 

flight in thermals with fairly fast glides between thermals. 

A portion of such a flight is sketched on Fig. 12. It 

includes a glide in a region of air descending at constant 

speed The sailplane flies at horizontal speed V, its 

rate of sink with respect to the surrounding air is v^, it 

covers the distance d in a time t^, losing height h. The 

sailplane then reaches a thermal assumed to have the shape 

of a vertical cylinder. The pilot begins circling flight, 

and the height h is recovered in the thermal in a time t^ . 

The rate of climb in the thermal is obtained by combining 

the sailplane rate of sink in circling flight v with the 

s o 

rate of climb of the thermal air v., which is assumed to 

th 

depend on the distance r between the sailplane and the 
thermal axis. 


V 


= v^j^(r)-v 


SC 


(A. 19) 


The cross-country speed is defined as the distance covered 
divided by the time required to reach the thermal and 
recover the altitude lost. 


^cc = 


(A. 20) 
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Figure 12: Model of flight used in the definition of cross- 

country speed. 
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= d/V = h/(v +v_, 
1 ' ' ' s ith 



(A. 21) 


so that 



(A. 22) 


The concept of cross-country speed was first introduced by 
K. Temmes (Ref. 30). It has been the starting point for a 
large number of investigations in the areas of sailplane 
design and flight strategy (see Refs. 31-33 for example). 

A. 3. 3 Flight Strategy 

Aside from keeping his sailplane as aerodynamically clean as 
possible, there are three major actions that a sailplane 
pilot may take to improve his cross-country speed for given 
atmospheric conditions. 

i) Achieve the best rate of climb in the thermal. 

Intuitively, the best value of will be obtained 

with the best rate of climb in thermals v . The 

c 

pilot must carefully select his position in the 
thermal to obtain the best possible lift. 

ii) Select the best interthermal speed. From the 


development on the speed polar (Eq. A.9), we know 


Ill 


that the rate of sink in straight flight is 

directly related to the horizontal flight speed V. 
For a given interthermal downdraft and a given 

anticipated achieved rate of climb, the pilot is 
left with the choice of the horizontal speed. 
Intuitively, if the next thermal is expected to be 
strong, the pilot should dive toward it at high 
speed. This may cause him to lose a lot of 
altitude in interthermal flight, but that loss will 
be recovered quickly in the strong thermal, and the 
cross-country speed will be high. Alternatively, 
if the pilot anticipates a weak thermal, he should 
choose a slow speed in interthermal flight, 
iii) Adjust the sailplane weight. To increase the 
weight of a sailplane alters its speed polar in 
such a way that a given rate of sink occurs at a 
higher horizontal flight speed. This is favorable 
in interthermal flight. However, the price to pay 
for that is a lower achievable rate of climb in a 
given thermal. It turns out that heavy weight is 
advantageous in typical strong thermal conditions 
but disadvantageous in typical weak conditions. 
Therefore, sailplanes are equipped with water tanks 
(water ballasts) in the wings. The pilot may 
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select the amount of water to carry on the basis of 
the thermal conditions for the flight. 

In the next two subsections, we will discuss how to 
evaluate the best achievable rate of climb for a sailplane 
in a given thermal and the best cross-country speed for 
fixed rate of climb and interthermal downdraft. 


A. 4 BEST RATE OF CLIMB IN THERMAL FLIGHT 
A. 4.1 Thermal Model 

A difficulty in any study of sailplane performance is the 
choice of the thermal model. We will follow Carmichael 
(Ref. 31) and assume the thermals to be cylinders with 
vertical axes. The distribution of vertical airspeed is 
given as a function of the distance r to the thermal axis by 

2/ ,r \ 

= cos I j if <2 

^thmax \ ^1/2/ ^1/2 


r 

= 0 if 52 (A. 23) 

^ 1/2 


where v., 

thmax 

is the 

maximum 

vertical 

airspeed 

( on the 

thermal 

axis ) 

and r ^/2 

is 

that 

radius 

where the 

vertical 

airspeed 

is 

reduced 

to 

one 

half 

of V. , 

thmax 

This 


distribution is plotted in Fig. 13, where typical values are 
given for the parameters. 
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Figure 13; Thermal models (arranged from Ref. 31). 


I 
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A. 4. 2 Best Rate of Climb 

From Eqs . A. 18 and A. 19 we obtain the sailplane rate of 
climb . 


For a given sailplane (fixed weight) in given atmospheric 
conditions, there are only two variables that the pilot can 
control to improve his rate of climb: the radius of turn r 
and the lift coefficient To evaluate the best 
achievable rate of climb, we specify 






(A. 24) 


[C^-(K/gr)2]^/^ 


c 


c 


0 


(A. 25) 


3r 


3C, 


L 


This yields the following set of equations 



0 


(A. 26. a) 


K 


r 


(A.26.b) 


^ +C_,C^ ) 


'DO L D1 L D2 L 
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These equations are nonlinear and must be solved 
numerically. 

Experiments with typical sailplane parameters have shown 
that the solution may yield an unrealistically high lift 
coefficient. When that happens, the lift coefficient must 
be set to its maximum allowable value C,. and Eq. A. 26. a 

solved for r. 

A. 5 BEST CROSS-COUNTRY SPEED 

The cross-country speed is given by Eqs. A. 9 and A. 22 
combined 


For a given sailplane (fixed weight) and known rate of climb 
and interthermal downdraft the pilot controls his cross- 
country speed by adjusting the horizontal flight speed. To 
evaluate the best available cross-country speed, we set 


Vv 


c 


V 


(A. 27) 


cc 



9V 


cc 


0 


(A. 28) 


3V 


This yields 
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C. 


DO 



V 


4 


V - C. 


0 


(A. 29) 


D2 


K' 


,2 


2K 


To obtain the best cross-country speed a given sailplane may 
reach in given atmospheric conditions, we proceed as 
follows ; 


use in the thermal from Eq. A. 26. 

ii) Find the resulting rate of climb from Eq. A. 24. 

iii) Determine the horizontal speed to select in 
interthermal flight from Eq. A. 29. 

iv) Find the resulting best achievable cross-country 
speed from Eq. A. 27. 


i) Determine the radius r and lift coefficient to 

Li 


Appendix B 

AERODYNAMIC SUBPROBLEM 

The objective of the aerodynamic subproblem is to design the 
wing so as to obtain the aerodynamic characteristics 
requested by the performance subproblem. These 
characteristics are the wing area, the coefficients of the 
drag polar, and the maximum sailplane lift coefficient. The 
design is subjected to the constraint that stall must begin 
at an inboard station of the wing so as to preserve lateral 
control in high angle of attack conditions. The variables 
of this subproblem are coefficients of the spanwise 
distributions of wing chord, thickness and geometric twist. 

The calculations of the contribution of the wing to the 
aerodynamic characteristics will be performed using accurate 
prediction methods since the wing design is the object of 
this subproblem. In particular, finite span and Reynolds 
number effects will be accounted for while estimating wing 
lift and drag (as suggested in 34). However, the wing will 
be assumed rigid, neglecting aeroelastic lift redistribution 
and similar phenomena. The effect of the remaining 
sailplane components on the aerodynamic characteristics of 
interest will be found using more approximate methods. This 
will enable us to keep the computational cost to a minimum 
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while still obtaining realistic estimates for these 
parameters . 

This appendix describes the calculations for obtaining 
the sailplane polar curve, its maximum lift coefficient and 
for estimating the position of the wing station where stall 
begins. In addition, it details the developments giving the 
control vectors A, H, and (see Chap. IV). 

B. 1 NOMENCLATURE 

{A} : vector of aerodynamic load parameters. 

{A^}, {A^}, {A^} : subvectors of A. 

Aj^ : horizontal tailplane aspect ratio. 

A^ : fuselage maximum cross-sectional area. 

[AS] : aerodynamic matrix, 

b : wing semi span. 

bj^ : spanwise position of break in wing planform. 

b^ : fuselage maximum width, 

b^ : width of wing element i. 

[B], [Blj, [B2] : Eqs. B.18 and B.20. 

c; c , c. , c. ; c. : section chord; chord for root, break 

r b t 1 

and tip sections; chord for wing element i. 
c :mean aerodynamic chord. 

c , , c , . , c , T : airfoil profile drag, minimum profile 
dp dpmin dpi 

drag, profile drag for unit lift coefficients. 

» » f 

^dpo' ^dpl' ^dp2 


: Eq. B.41. 
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c, , c, . , c. / c, : airfoil lift coefficient, design lift 
1 li Imax la 

coefficient, maximum lift coefficient, lift curve 
slope . 

I » f 

c, , c, , c, : airfoil three-dimensional lift, lift 

1 lo la 

for zero fuselage angle of attack, lift curve slope 
coefficients . 

c : airfoil coefficient of moment about aerodynamic 

mac 

center . 

C C ^ . : coefficients of aerodynamic moment for wing 

mac macw 

fuselage combination, for wing alone. 

: coefficient of aerodynamic forces acting on 

sailplane parallel to mean aerodynamic chord. 

^D' ^Dw' coefficients of drag for sailplane, for wing. 

C^S : drag area. 

Cp : skin friction coefficient. 

^L' ^Lo' ^La ■ trimmed lift, trimmed lift at zero 

fuselage angle of attack, untrimmed lift curve 
slope . 

^Lh' ^Lha " horizontal tailplane lift, lift curve slope 
coefficients . 

*^Lw' ^Lwo' *^Lwo ' lift, zero fuselage angle of attack 

lift, lift curve slope coefficients. 

^Lwf ■ coefficient for wing fuselage combination. 
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*^N' ^Nw " aerodynamic forces acting normal 

to mean aerodynamic chord on sailplane, on wing, 
dj^, dj^ : chordwise and normal moment arms of aerodynamic 

forces with respect to axis Y (Fig. 18). 

^aN' ^aC ' forces per unit length acting normal 

and parallel to mean aerodynamic chord (Fig. 18). 
hj^, h^ : Fig. 16. 

H, : vectors describing shape of wing, shape of wing 

element j . 

i^ : wing angle of incidence (Fig. 14). 


Kj, Kjj: 

^f' ^fn' 


m 


aY 


n 


n. 


n 


n. 


Hxt / n. , ; 

N Nw 
Re : 

S : 


Eq . B . 26 . 

Ij^ : Fig. 16. 

aerodynamic moment per unit length about spanwise 
axis . 

number of wing stations in integration of lifting 
line equation. 

longitudinal sailplane load factor. 

number of wing elements at local structural level, 
number of intervals in discretization of spanwise 
distributions of aerodynamic load parameters, 
normal sailplane, wing load factors. 

Reynolds number, 
wing planform area. 
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S^, ^fwet ' fi^selage planform area, wetted area, 

portion of wing area inside fuselage outline (Fig. 

16 ). 


S^, : planform areas for horizontal and vertical tail, 

t; t.^, tj^, t^; t^ : wing thickness; thickness for root, 

break and tip sections; thickness for wing element 
i . 


: fuselage volume. 

; indicated airspeed. 

W : sailplane gross weight. 

^ac' ^cg' ^n *' ^oi^gi'tudinal position of aerodynamic center 
of sailplane without tail, of sailplane center of 
gravity and neutral point; datum is at leading edge 
of mean aerodynamic chord. 

X , X . : dimensional chordwise position of airfoil 

ac eg mt ^ 

aerodynamic center, area center of gravity, maximum 

thickness; datum is at leading edge. 

jX } : vector of variables of aerodynamic subproblem. 

; vector of performance variables. 

y : dimensional spanwise coordinate. 

y : spanwise position of mean aerodynamic chord. 

z , z , z . : dimensional position of airfoil aerodynamic 

ac eg mt 

center, area center of gravity, maximum thickness 

normal to airfoil chord; datum is at airfoil chord, 
positive up. 
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p 

R 

O 

“or 

a , , a 

olr 

“ol^at 

: Eqs . 

B.24. 



section 


C/ Tl/ Y: z 


section geometric twist; geometric twist for 
break and tip wing sections. 

geometric twist at section of mean aerodynamic 
chord. 

normalized with respect to the 


Tl 

^h 

P 

0 

X, 6 


3: 


ac' eg' mt 
section chord. 

non-dimensional spanwise coordinate (y/b). 
horizontal tailplane taper ratio, 
mass density of air. 

non-dimensional spanwise coordinate (Eq. B.15). 

normalized with respect to section 


ac 

chord 


V : kinematic viscosity of air. 

Commonly used superscripts 

1, 2, 3 : aerodynamic load parameters related to the 

various structural design requirements. 

' : three-dimensional airfoil aerodynamic quantity 

(i.e., a quantity corrected for finite span 


effects ) . 

T : transposed vector or matrix. 

-1 : inverse matrix. 

Commonly used subscripts 
i as in : element i of vector A, or 
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: discretized value of function A(y), i.e., = 

A(Yi). 

ij : element ij of matrix, 

b : wing section at planform break, 

f : fuselage, 

h : horizontal tail, 

m : miscellaneous, 

r : wing section at root, 

t : wing section at tip. 

V : vertical tail. 

wf : wing fuselage combination. 

Special symbols 
{ } as in {A} : vector A, or 

: discrete values of function A(y) collected 
: in vector {A}"^ = { A( y^ ) , A( y^ ) , . . . j . 

[ ] : square matrix. 


B.2 WING DESCRIPTION 
B.2.1 Wing Shape 

The wing aerodynamic model is described in Fig. 14. It is a 
straight cantilever wing with linearly varying chord c, 
thickness t and geometric twist t. The wing is positioned 
so that the chord of the root section is at an angle i 
(wing angle of incidence) with respect to the fuselage axis. 
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The chord of the wing section at station y makes an angle 
e(y) with respect to the chord of the root section. The 
sections are arranged so that the line joining their point 
of maximum thickness is straight and perpendicular to the 
sailplane plane of symmetry. 

In terms of the spanwise coordinate y, we have 
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Zero-lift line 
» section y 



Figure 14: Wing definition. 
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The quantities b, h^, c^, c^, t^,t^,t^, and e^ are 
the design variables of the aerodynamic optimization 
subproblem; they are collected in vector X^. 

The total wing area is 


S = b, (c +c, ) + (b-b, ) ( c, +c . ) 
b'r b^' b^'b t' 


(B.4) 


The wing mean aerodynamic chord is defined by 

2 2 , . , 

c = -g / c (y)dy 


(B.5) 


The spanwise position of the aerodynamic center is obtained 
from 

'b 


y = 4- / c ( y ) ydy 


(B.6) 


while the wing twist at the aerodynamic center is found from 
Eq. B.3 as 


t = E(y) 


(B.7) 


B.2.2 Airfoil Description 

In the past thirty years, sailplane wings have been built 
with specially tailored, laminar-flow airfoils (see Ref. 
35). However, for this application, we will use airfoils 
from the NACA 6-series, a class typical of immediate post 
World War II sailplane design. This choice was conditioned 
by the availability of experimental data. The specific 
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airfoil family chosen is the NACA 64-4XX which is 
extensively described in Ref. 36. For the purpose of this 
application, linear regressions were developed from the data 
collected in Ref. 36. These laws relate the airfoil 
performance characteristics to their thickness ratio (t/c) 
and the Reynolds number based on the chord (Re=Vc/v). The 
data from Ref. 36 is valid for the following ranges of 
thickness ratio and Reynolds number 

.09<t/c<.18 and 630000<Re<6300000 (B.8) 

Typical Reynolds number for sailplanes may be as low as 
400000; therefore, the laws are somewhat deficient. 
However, no sufficient set of data was found for Reynolds 
numbers less than 630000. 

The following performance laws are used (see Fig. 15). 
i) Lift coefficient (c^) versus angle of attack (o^)- 


= (-98. 13+12. 4291ogRe)+(1964. 96-249. 9211ogRe) (t/c) 



c, <c, 

1 Imax 


0 


^'l^^lmax 


+ (-13785. 2+1722. 21ogRe) (t/c) 


2 


+( 32757. -4056. logRe) (t/c) 


3 


c 


la 


= (- . 12661+.035091ogRe)+(1.1562-.165711ogRe) (t/c) 
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a) Lift coefficient. 




c) Moment Coefficient 
(about aerodynamic center). 


d) Definition of geometry 
and aerodynamic coefficients. 


Figure 15; Wing section, definition of geometry and 
aerodynamic performance coefficients. 
Each curve corresponds to one Re and 
t/c combination. 
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*^lmax “ ( 1. 756-. 20581ogRe) + (-66. 826+12. 8401ogRe ) (t/c) 
+ (306. 11-56. 7461ogRe ) (t/c)^ 


(B.9) 


ii) Moment coefficient about aerodynamic center (c ) 

mac 


c = -.075+.0556(t/c) 

mac ' ' ' 


(B.IO) 


iii) Profile drag coefficient (c^^) versus lift coefficient. 


^r^ii 


'"dp ^dpmin"*^ ^ ''dpi ^dpmin^ 


1-c 


li 


Cdpmin ~ ( •1315-.0411351ogRe+.00328131og Re) 


+( . 1956- .0272861ogRe) (t/c) 


'dpi 


= (8.026-2.2487961ogRe+. 15737041og Re) 


+( -126. 0558+33. 683971ogRe-2 .228871og Re) (t/c) 


+(583 .6108- 139. 02771ogRe+7. 8373 log^Re) (t/c) ^ 


+ (-656. 255+96. 2961ogRe) (t/c)' 


^li 


.4 


(B.ll) 
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iv) Aerodynamic center position (x and z ). 

3.C 3C 


ac 


= xc 


-6 . 2579+2 . 79581ogRe- . 42021og^Re+ . 0211231og^Re 


+(9.308-.4941ogRe+.02231og Re) (t/c) 


+(-54.81+.281ogRe) (t/c) ^+130.7 (t/c)^ 


z = Cc 

ac 


c = -1.8l43-.19111ogRe+.16271og^Re-.0147571og^Re 


+( 69. 045-18. 325 logRe+1.33781og^Re) (t/c) 

+(-94. 67+5.241ogRe) (t/c)^+164.6(t/c)^ (B. 12) 


v) Maximum thickness point position (x ^,z ^). 

^ ^ ' mt' mt' 


^mt 


pc p = .37 


^mt 

= 

Zc y = .02 

(B. 13 ) 

vi ) 

Area 

center of gravity position • 


^cg 

= 

6c 6 = .42 


z 

r* rr 

= 

nc Ti = . 02 

(B.14) 
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B.3 LIFT AND AERODYNAMIC MOMENT CALCULATION 
B.3.1 Wing Lift 

The three-dimensional spanwise distribution of lift 
coefficient may be obtained by integration of Prandtl's 
lifting line equation. This classical integral equation is 
developed in Refs. 37 and 38 where it is said to be 
applicable to straight, high aspect ratio wings in 
imcompressible flow. Several solution schemes exist for 
Prandtl's equation. We will use Glauert's method, which is 
based on Fourier series. 

A new coordinate is defined along the wing by 

(t> = cos"^(y/b) (B. 15) 


where the angle <t> varies between 0 at the wing tip and tr/2 
at the wing root. Stations are chosen along the semispan at 
equal interval L<t> . If there are n intervals 

Si 


A0 = ir/(2n ) 

a 


(B.16) 


The three-dimensional lift coefficient (c^^) at an 
arbitrary section along the wing is given by 


' 2b 

C-, = 


n 




1 = F(J) B sin[(2j-l)«] 
' ' 3 = 1 


(B.17) 
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For a symmetric lift distribution, 
obtained from 


the coefficients B . are 

J 


{B1 = [AS] [a] 


[AS] 


-1 

ij 


2b 


c, (iA0)c(iA^) 
L_ lot 


(2j-l) 


4sin( iA<^ ) 


sin[ i (2 j -1 ) A0 ] 


= a(iA(^) 


(B.18) 


a(iA0), the angle of attack for section i, is measured with 
respect to the zero lift line of that section. The lift 
curve slope of the airfoil c^^{iL(t>) is a known function of 
the section chord c(iA0) and thickness t(iA0) defined in 
Subsec. B.2.1 and, also, of the Reynolds number 
(Re=Vc(iA^)/v) . 

From Fig. 14 we see that 


a. = a,-+o. 
1 f 1 


a. = i +E.-0 

1 w 1 oil 


If we define 


(B.19) 


[Blj = 

[AS] [u] 


[B2] = 

[AS] [a' ] 


T 

[u]^ = 

[1,1 1} 

(B.20) 
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we obtain the three-dimensional wing lift distribution as a 
function of the sailplane angle of attack 


Ct ( y) = c, ( y )a -+c, 

1 \ JT / la' ' f lo 


n 


“ cfyj j^l BljSin[ (2j-l)0(y) ] 
n 

^ cf^ j=l B2^sin[ (2j-l)0(y) ] 


(B.21) 


where, of course, 0(y) is given by Eq. B.15. 

The total wing lift coefficient is now obtained from 

2 

^Lw = tI 

or, using Eq. B.21 


C, = Ot,r+C^ 

Lw Lwa f Lwo 


"^LWa = Bl^)/S 


""lwo = B2^)/S 


(B.22) 


For subsequent calculations, it is convenient to write 
the wing lift coefficient as 


C_ = C. (a_ + i -a , -a -e . ) 
Lw Lwa f w olr olat' 


(B.23) 


Comparing Eqs. B.22 and B.23, we obtain 

C, 


ol at w olr C 


Lwo 


Lwa 


(B.24) 
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The quantity is called the change in zero-lift angle of 

attack of the wing per unit of positive aerodynamic twist at 
the wing tip, being the aerodynamic twist at the wing 

tip. It must be emphasized that in order to complete these 
calculations, the Reynolds number Re and the fuselage angle 
of attack must be known. This point will clarified in 
Subsec. B.4.3 

B.3.2 Sailplane Trimmed Lift 

In Ref. 39, the total trimmed lift is given as, 

/ X -X 

/ eg 

= ^Lwf P 

\ Ih 

The lift coefficient for the wing fuselage combination 

the distance between the sailplane center of gravity and its 

aerodynamic center aerodynamic moment 

coefficient about the sailplane aerodynamic center for the 

complete aircraft minus horizontal tail c will be 

mac 

calculated in this section. The mean aerodynamic chord c 

was defined in Subsec. B.2.1, and the tail length 1, is a 

n 

fixed parameter of the design problem. 

The lift coefficient for the wing-fuselage combination is 


ac 


mac 


(B.25) 


given for a mid-wing configuration by 
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'Lwf 


K 


II 


[a^-a tE .+ 

I Lwa f ol at 


i -a , ) ] 

w olr' 


K, 


2-irb^ 

Kj = (l+2.15b^/b)(l-S^yS) g 

Lwa 


Kjj = (l+.7b^/b) (1-S^yS) (B.26) 

The quantities b^ and are fixed parameters in this 
design problem and they are defined in Fig. 16, The 
longitudinal distance between the aircraft center of gravity 
and its aerodynamic center (tail effects not included) may 
be written 


X -X =(x -X )+(x -X ) 
eg ac ' eg n' ' n ac ' 


(B.27) 


where x^ defines the stick-fixed neutral point, i.e., that 

position of the center of gravity for which the sailplane is 

neutrally stable for constant elevator angle. The quantity 

(x -X )/c is called the static margin; it is a measure of 
n eg 

longitudinal static stability and should be positive for 
stability. It may be taken as a fixed parameter for the 
problem. The quantity obtained as follows 


for an aircraft with a T-tail 
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Figure 16: Geometric description of the sailplane. 
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n ac 


c 


^Lha 


C 


La 
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C 


La 


+C_ 

I Lwa Lha 



(B.28) 


The tailplane lift curve slope may be estimated by 

2it 


'Lha 


2(l+2Xj^) 


(B.29) 


1 + 


A^(l.Xh) 


where Aj^, X^^ and are the tailplane aspect ratio, taper 
ratio and total planform area; they are fixed parameters. 
The downwash gradient de/da is found from 



(B.30) 


where b is the wing semi span, Ij^ and hj^ are dimensions 
related to the tail (Fig. 16). To trim the sailplane, the 
horizontal tailplane must exert a lift given by the 
coefficient 


'Lh 


Sc 


^h^h 


c +C^ J 

mac Lwf 

/X -X \ 

f eg ac 




1 

i ^ 1 



(B.31) 
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B.3.3 Pit-ching Moment Coefficient 

For the sailplane without horizontal tail, the pitching 
moment coefficient is given by 


C = C +A-C 
mac macw f mac 


(B.32) 


The wing moment coefficient is given for straight wings by 
'b 


macw 


where c 



c (y)c (y)dy 
mac ^ ^ 


mac 


(B.33) 

is the airfoil pitching moment coefficient. It 


is a known function of the thickness ratio (t/c) as 
explained in Sec. B.2.2. 

The effect of the fuselage on the pitching moment 
coefficient is (see Ref. 39) 


b, nb^h^-l-CT 

A.C = - .9(1-5. 

f mac 1-' ScKt-C^ 

f I Lwa 


(B.34) 


All the elements in Eq. B.34 were calculated before or are 
defined in Fig. 16 except for the total sailplane lift 
coefficient at zero fuselage angle of attack («£)• 
Combining Eqs. B. 25-26, B.32 and B.34, we obtain 

K, 


KtC, 

I Lwa 


^ 'II / • X 

(1 “01 1 ) 

ol at Kj ' w olr' 


-X 




1, macw 
h 
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B.3.4 Maximum Lift Coefficient 

The information obtained in the previous sections may be 
used to estimate approximately the sailplane maximum lift 
coefficient . 

The three-dimensional distribution of wing lift 
coefficient is given by Eq. B.21. Assuming that stall 
occurs in a given wing section when that section lift 
coefficient equals its maximum local lift coefficient 
obtained from experimental data (Eq. B.9), we may define the 
fuselage angle of attack for which stall occurs in section y 
as 


The wing section where stall initiates (at spanwise station 

y^) is such that o£g(y) is minimum. It may be determined by 

numerical minimization of “fg(y) with respect to y. The 

resulting angle of attack may be used with Eq. B.25 to 

determine the maximum sailplane lift coefficient 

^ Lmax 

t 

This approach is described in Ref. 39. The 

f 

Cio(y) and ^]_niax^^^ distributions are Reynolds number 
dependent, as was pointed out earlier. Therefore, the 


procedure should be repeated iteratively. Such a refinement 
does not seem appropriate in view of the approximate 
character of the method; the above mentioned distributions 
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will be calculated with Reynolds numbers based on a typical 
(slow) flight speed. 


B. 4 SAILPLANE DRAG PREDICTION 
B.4.1 Wing Drag 

In this subsection, we develop expressions for the vortex 
induced drag distribution, the total drag distribution and 
the total wing drag. 

The local vortex induced drag coefficient is given in 
Ref. 38 as 

Cdi(y) = c^(y) [a(y)-c^(y)/c^^(y) ] (B.37) 

I 

where the ( ) emphasizes the fact that the vortex induced 

f ! 

drag coefficient c^^ and the lift coefficient c.^ include 
finite span effects. Now, from Fig. 14 


“(y) = af+i^+E(y)-“oi(y) 


(B.38) 


The total drag distribution on the wing is 


c^(Y) = 


(B.39) 


where the local airfoil profile drag coefficient c^^ is 

f 

known once the local lift coefficient c^, the local 
Reynolds number and the local thickness ratio are known, as 
shown in Subsec. B.2.2. 


The total drag on the wing is found as 
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\ 


b , 



(B.40) 


It may be calculated once the flight speed, the fuselage 
angle of attack and the wing geometry are known. 

For the calculation of the loads on the wing, it will 
prove necessary to know the wing drag as a function of the 
fuselage angle of attack. Using Eqs. B.ll, B.21, and 
B. 37-40, we obtain 








+c, (y) 

dpmin ' ^ ' 



2 





1 . -c 


li 











2 
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l.-c 


li 


(B.41) 


B.4.2 Total Sailplane Drag 

The various contributions to the sailplane drag from sources 
other than the wing will be given in this subsection, 
following Ref. 39. 

B.4.2.1 Fuselage Drag 

The drag of the fuselage includes viscous and induced 
components . 



F fwet 



2.2 3.8 

+ 


eff ^eff 


1 


f 


X 


eff 


(B.42) 


(4A^A)^/2 



the fuselage volume, wetted area, planform area, maximum 
cross-sectional area and length; they are fixed parameters 


143 


ORJGSfv.AL P/ttSS 
OF POOR QUALITY 

of the design problem. The angle is the fuselage angle 
of attack. is the skin friction coefficient; it is a 

function of the Reynolds number and the location of the 
boundary layer transition on the body considered. For lack 
of sufficient information, we will assume that transition 
occurs at the nose of all sailplane components, a 
conservative assumption. Then 

.4550 

Cp = (B.43) 

, , . 2 . 58 

( logRe ) 

For the fuselage, the Reynolds number is based on the 
total fuselage length (l^). 

B.4.2.2 Wing-Fuselage Interaction Drag 

The drag due to wing- fuselage interaction may be obtained 
for straight wings from: 

.138(bp/b)C^^S^ 

^V^wf = (2-nbp/b) + 6.75CpC^t^ (B.44) 

(l-Hc^/c^)Trb^ 

In this equation, c^, t^, b, and S are, respectively: the 

wing root chord, root thickness, semispan and total area. 
The lift coefficient for zero fuselage angle of attack 
was obtained in Subsec. B.3.3. Also, b^ is the fuselage 
radius. The skin friction coefficient is obtained from 
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Eq. B . 43 , with the Reynolds number based on the nose-wing 
distance . 


B.4.2.3 Horizontal Tailplane Drag 

Combining the drag due to the tailplane lift and the 
tailplane viscous drag, we find 


(C^S)j^ = 1.35 


s 

^Lh^h 


irA, 


+ 2Cp[l+2.75(t/c)j^]S^ 


(B.45) 


^h' ^h' respectively: the tailplane planform 
area, aspect ratio, and thickness ratio; they are fixed 
parameters in this problem. The tailplane lift coefficient 
Chh was obtained in Subsec. B.3.2. The skin friction 
coefficient and the thickness ratio are based on the average 
tailplane chord. 


B.4.2.4 Vertical Tailplane Drag 

The main contribution to the vertical tailplane drag is 
viscous in nature, it is given by; 

(CpS)^ = 2Cp[l+2.75(t/c)^]S^ (B.46) 

and (t/c)^ are the vertical tailplane planform area and 
thickness ratio. The skin friction coefficient as well as 
the thickness ratio are based on the average vertical 
tailplane chord. 
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B.4.2.5 Miscellaneous Drag Contributions 

A sailplane is aerodynamically clean in design. The only 
drag contribution that will be considered in addition to 
those already mentioned is that due to the gaps of the flaps 
and control surfaces. For movable surfaces along the entire 
trailing edge of the wing and the vertical tailplane, we 
have : 

(Cj^S)^ = 2(S+S^)10'^ (B.47) 


B.4.3 Drag Polar Calculation 

The total sailplane drag coefficient is now given by 


S =s'<Vw"<V>£"'S^'w£*<SS)h*<SS>v*<Vm> 


The various contributions were obtained in the two previous 
subsections . 

To evaluate the drag polar (Subsec. A. 2.1) for known 
weight and geometry, we proceed as follows. 

i) Choose a lift coefficient value and determine an 
approximate flight speed by equating the lift to 
the known weight (Eq. A. 8). At this point, the 
local Reynolds number is known for all sailplane 


components . 
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ii) Proceed with the calculations of Subsec. B.3.1 and 

obtain all the information on the wing lift. The 
fuselage angle of attack is still unknown at 

this step. 

iii) Calculate all the terms in the sailplane trimmed 

lift equation (B.25), as shown in Subsecs. B.3.2-3. 

iv) Find the fuselage angle of attack from the 

information obtained in iii) and the lift 

coefficient selected in i). 

v) Continue with the calculations of the various drag 
components as described in Subsecs. B.4.1-2. 

The drag polar is obtained point by point by repeating 
this procedure for different lift coefficients. In this 
application, the curve is modeled by a parabola; therefore, 
the process must be repeated only three times to obtain the 
coefficients in Eq. A.l. 

B.5 CONTROLS CALCULATION 
B.5.1 Wing Shape (Vector H) 

In order to design the structure of the wing, the global 
structural subproblem must have access to the wing shape. 
The shape information to send from the aerodynamic 
subproblem to the global structural subproblem is collected 
in vector H 

{Hj^ = {b,b^,Cj.,Cj^,c^} 


(B.49) 
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To obtain a linear approximation to the constraint 
satisfactions in the global structural subproblem, we need 
the matrix 3H/8X (see Chap. IV). This matrix is easily 
obtained from Eq. B.49. 


When performing the detailed structural design, the wing is 

modelled by n cylindrical elements of equal length and 
s 

constant geometric and structural properties. The element 
properties are chosen to be those of their middle section in 
the actual wing. Three shape parameters must be known for 
the analysis of any one element. They are given for element 
j by the vector 


The reference station for wing element j is determined by 


B.5.2 Wing Element Shape (Vector ) 



T 



(B.50) 


where b ^ , c^and t^ are, respectively the element length 
chord and thickness. The element length is 


b . = b/n 
3 ^ s 


(B.51) 


If we define 


m^ = (2j-l)/(2ng) 


(B. 52 ) 



(B.53) 
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The element chord (c.) and thickness (t.) are now found from 

Eqs. B.1-2, where y is set to . 

While optimizing the aerodynamic subproblem, we need 

linear approximations to the constraint satisfaction in the 

local structural subproblems in terms of the aerodynamic 

variables (see Chap. IV). These linear expansions contain 

the matrix 3H./3X , which may be obtained from Eq. B.50. 
J « 

Its non zero entries are as follows. 


3bj/3b 

3Cj/3b 

3c ./3b. 

1 b 

3c ./3c 
y r 

3c ./3c, 

1 b 


^^3 

[m^(Cb-Cr)]/bb 
[ (l-m^)(c^-Cj^)bj^]/(b-bj^)2 
[m^ (Cj.-Cj^)b]/bj^^ 

[ (1-m^ ) (Cj^-c^)b]/(b-bj^)^ 

3tj/3t^ = 1-m^b/bj^ 

0 

3tj/3tb = n>^b/bj^ 
l-(m^b-bj^)/(b-bj^) 


y 

t> 

^1 b 

y -^b, 
^1 b 

y • <b. 

^1 b 

y -^b, 
^3 b 
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3c ./3c. = 3t ./3t^ = 0 

I t 1'^ t 

y • <t>. 

^1 b 


= (m^b-bj^)/(b-bj^) 

^3 t) 


3t^/3b = [m^ (t^-t^) ]/bj^ 

y • <ti, 

^3 b 


= [ (l-m^)(t^-tj^)bj^]/(b-b^)^ 

^3 b 


3tj/3bj^ = Im^ (t^-tj^)b]/b^ 

y • <b. 

^3 b 


= ( (l-m^)(tj^-t^)b]/(b-bj^)^ 

y -^b, 
^3 b 

(B.54) 

B.5.3 Aerodynamic Loads Calculation 

The Federal Aviation Agency requires that the 

integrity of the sailplane be verified for a 

structural 
number of 


different load cases (see Ref. 40). The cases considered in 
this study correspond to symmetric wing loadings, and their 
selection will be discussed in App. C. Each symmetric load 
case is uniquely defined by an indicated airspeed and a 
wing load factor The total load on the wing is then 
given in coefficient form by 


n., W 
Nw 


Cm = 2 
Nw 


psv: 


(B.55) 


The total wing load is assumed to act in a direction 
perpendicular to the reference chord of the section 
corresponding to the mean aerodynamic chord, as shown on 
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Definition of force coefficients for the wing 
alone and the total sailplane. 


Figure 17 ; 
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Fig. 17. 

To determine the spanwise distributions of aerodynamic 
forces and moment we must find the fuselage angle of attack 
which corresponds to Following the developments 

outlined in Subsecs. B.3.1 and B.4.1, we obtain the wing 
lift and drag coefficients (C. and C_ ) as functions of 
(Eqs. B.22 and B.41), the Reynolds number being based on V^. 
Then it is clear from Fig. 17 that 


'Nw 


= (a^)cos(a^+i +e)+C^ (a^)sin(a-+i +e ) 
Lw' f' ' f w ' Dw' f' ' f w ' 


(B.56) 


In what follows we will assume that the load cases are 

selected so that a^+i +t is small, so 

r w 


'Nw 


= ^Lw(“f) 


(B.57) 


The only unknown in Eq. B.57 is the angle of attack a^, 
which may be easily solved for. Then, the spanwise 

t I 

distributions of lift c^^ and drag c^ are found from Eqs. 

B.21 and B.41, respectively. The distribution of 
coefficient of pitching moment about the section aerodynamic 
center is taken from Eq. B.IO, ignoring finite span 

effects. This situation is summarized on Fig. 18 for an 
arbitrary wing section. For subsequent calculations, we 
express the aerodynamic forces and moments in a system of 
axes centered in the plane of symmetry of the sailplane 
with: 


t 

I 
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4 Axis N 



Figure 18: 


Wing section: aerodynamic forces applied and 
corresponding geometric parameters. 
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i) axis Y previously defined, spanwise and 
perpendicular to the plane of symmetry, 

ii) axis C in the plane of symmetry and parallel to the 
mean aerodynamic chord, 

iii) axis N in the plane of symmetry also and 
perpendicular to axes Y and C. 

Note that axis Y joins the points of maximum thickness for 
each section. Those points are located at distances &c from 
the leading edge of the section and Zc above the chord line; 
3 and y are constants characteristic of the selected airfoil 
and are given by Eq. B.13. 

At section y, the forces per unit length reduce to (see 
Fig. 18) 

^aN^^^ = -5pV^c(y)c^(y) 

fac^y^ = -5pV^c(y)c^(y) (B.58) 

These forces act at the aerodynamic center of the section, 
which is located at distances x(y)c(y) from the leading edge 
and C(y)c(y) above the chord line. x and c are also given 
in Sec. B.l; like the other airfoil parameters, they are 
functions of the thickness ratio and the Reynolds number. 
Measured positively along the axes C and N, the moment arms 
of the forces with respect to axis Y are (see Fig. 18) 
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dc(y) = -[B-x(y) ]c(y) 


dj^(y) = -[3f-c(y) ]c(y) 


(B.59) 


The total moment per unit length of the aerodynamic 
forces at section y (positive clockwise) is finally 


We conclude this subsection by finding the total 

sailplane load factors. At this point we know the angle of 

attack yielding the wing load factor n., and, also, the 

flight speed and, therefore, the Reynolds number. We 

perform the calculations described in Subsecs. B.3.1-3 and 

B.4.1-2 to find the lift and drag coefficients C_ and C_. 

Lj d 

Then, from Fig. 18 the aerodynamic force coefficients in the 
YCN set of axes are (again for small (o^+i^+e)) 




C. 


N 


C 


L 


C 


C 


C. 


D 


(B. 61) 


The total sailplane load factors are 





2W 



n, 


C 


2W 


(B. 62) 
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B.5.4 Load Parameters (Vector A) 

The wing will be designed so as to meet three structural 
requirements (see App. C). The aerodynamic parameters 
needed to perform the structural design will be calculated 
in the aerodynamic subproblems and communicated to the 
global structural subproblem through vector A, which will be 
partitioned into three vectors, one per structural 
requirement. 


Requirement 1 limits the displacement of the wing tip. 
This necessitates the knowledge of the sailplane normal load 
factor n^^ and of the spanwise distribution of aerodynamic 
force normal to the mean aerodynamic chord The latter 
will be transmitted to the global structural subproblem as a 
vector containing the values of calculated at 
(n^+1) points along the wing; those points being equally 
spaced and including the wing root and tip. Therefore 



(B. 63 ) 



(B.64) 


Requirement 2 is that no failure or buckling occur in the 
wing elements in a high-speed positive pull-up maneuver. To 
perform the necessary calculations, the sailplane load 
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factors and as well as the corresponding spanwise 

distributions of aerodynamic transverse force f^j^ and moment 

m „ must be known. Hence 
aY 


,.2,T_ , 2 2 . J2. ,T , 2 

{A } f n^^ , n^ , I f , { m^Y 1 1 


(B.65) 


2 2 

where ^^aY^ obtained by the process 

described above. While structural design to requirement 2 

is the object of the local structural subproblems (App. D), 
2 

vector fA } must first be transmitted to the global 
structural subproblem so as to account for the effect of 
inertia on the spanwise distributions of forces and moment. 

Finally, requirement 3 is that the wing torsional 
divergence speed must exceed a given limit. As will be 
shown in App. C, the spanwise distributions of lift curve 

t 

slope c^^ and chordwise position of the aerodynamic 
center x are required for the calculation of the divergence 
speed. Therefore 


3 T ' 3 T 3 T 

\Pr\ = 


(B.66) 


' 3 3 

where, again, and {x } contain appropriate values 

f 

of c^^(y) and x(y) for equally spaced wing stations. 

The matrices 3A/3X and 8A/3X are needed when 

a P 

constructing linear extrapolations to the measures of 
constraint satisfaction in the global structural subproblem 
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(App. C). These matrices will be found by 
differences . 


finite 
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Appendix C 

GLOBAL STRUCTURAL SUBPROBLEM 

The global structural subproblem objective is to design the 
wing structure so as to achieve the total sailplane weight 
specified in the performance subproblem. Once a weight is 
chosen for the sailplane fuselage, tailplane, equipment and 
pilot, the available weight for structural and non- 
structural wing elements is known. The design variables of 
this subproblem are the coefficients in the spanwise 
distributions of weight and flexural and torsional 
stiffnesses. The optimization is subjected to the 
constraints that the wing tip deflection is limited in high 
speed straight flight and that the torsional divergence 
speed is beyond a set limit. 

The wing shape and aerodynamic load parameters are 
obtained from the aerodynamic subproblem. The wing flexural 
and torsional responses are investigated using Rayleigh-Ri tz 
analyses. The redistribution of aerodynamic forces due to 
the wing elastic deformations is neglected. 

Optimization and sensitivity analysis require that 
derivatives of the wing weight, tip displacement and 
divergence dynamic pressure be calculated with respect to 
the subproblem variables and parameters. Unless otherwise 
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specified, these derivatives are found by finite 
differences . 


C. 1 NOMENCLATURE 

[A| ; vector of aerodynamic load parameters. 

{A^}, {A^j, {A^} : subvectors of |A}. 

b : wing semi span, 

c : wing chord. 

I f 

c, , c, : section three-dimensional elastic lift, lift 

le la 

curve slope coefficients, 
e^, * Fig. 19. 


El, EI^, : distribution of wing bending stiffness, ith 

coefficient in that distribution (Eq. C.l), target 
bending stiffness for element j . 

EI^, EI^ : bending stiffness distribution for the linear 
wing model, root bending stiffness for that model. 

I 

^aN' ^aNi' ^aN * aerodynamic forces normal 

to the mean aerodynamic chord, aerodynamic force at 
spanwise station i, distribution estimated by 
Lagrangian interpolation (Eq. C.4). 

f., : distribution of total forces normal to the mean 

N 

aerodynamic chord. 

f^, f^ : distribution of transverse forces on the linear 

model, root value of that distribution. 
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: vector containing the end forces acting on wing 

element j . 

{H] : vector containing parameters describing the wing 

shape . 


GJ, G J . , GJ . : distribution of wing torsional stiffness, 

1 0 J 

ith coefficient in that distribution (Eq. C.l), 
target torsional stiffness for element j . 

Kj : vector containing the target stiffnesses for 

element j . 

m : distribution of moment of forces about spanwise 

axis including effect of elastic deformations. 


m 


aY' 


m 


eaY 


m 


aYi 


, m^Y = distribution of aerodynamic moment about 
spanwise axis, aerodynamic moment at spanwise 
station i, distribution estimated by Lagrangian 
interpolation (Eq. C.4). 

distribution of elastic aerodynamic moment about 
spanwise axis. 


* * 

m m : Eqs. C. 16-17. 

eaY eaYi ^ 

m^ : distribution of total moment of forces about 

spanwise axis, includes aerodynamic and inertia 
forces . 

: bending moment in section i . 
n^^ : sailplane longitudinal load factor. 
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number of intervals in the discretization of the 
distributions of aerodynamic characteristics, 
sailplane normal load factor, wing normal load 
factor . 

number of intervals in the discretization of the 
distributions of wing weight and stiffnesses at the 
global level, number of wing elements in the 
discretization of the wing at the local level. 

: Lagrange polynomials, Eqs. C.l, C.4. 
dynamic pressure, divergence dynamic pressure, 
ultimate safety factor, 
torque at wing station i. 

^tmax ■ bending mode, non-dimensional 

bending mode, bending deflection at the wing tip, 
maximum value for u^. 

Vj^, ^Dmin ' indicated flight speed, maximum 
indicated dive speed, divergence speed, minimum 
divergence speed, 
shear force at wing station i. 

w^j : spanwise distribution of wing weight per unit 
length, ith coefficient in that distribution (Eq. 
C.l), target weight per unit length for element j. 
sailplane gross weight, 
weight of sailplane without wings. 
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y, Yi 

N 

: spanwise coordinate, spanwise coordinate for ! 

station i . 

{X 1 
* g 

: vector containing global structural subproblem 

design variables. 

a 

: slope of the force distribution on the linear wing 

model ( Eq . C . 9 ) . 

a 

r 

: wing section rigid angle of attack. 

e 

: slope of the stiffness distribution on the linear 

wing model (Eq. C.9) 

B, 6 

: chordwise position of the section elastic axis, and 
center of gravity with respect to section leading 
edge, normalized with respect to section chord 
(Fig. 19). 

E , E 

: wing geometric twist (see Subsec. B.2.1), wing 

twist at the mean aerodynamic chord. 

T1 

: normal position of the section elastic axis and 

* 

T1 , T| 

center of gravity with respect to section chord, 
normalized with respect to section chord (Fig. 19). 
: wing bending modes. 

X; 

: chordwise position of section aerodynamic center 

with respect to the' section leading edge, 

normalized with respect to the section chord; value 
of X at station i. 

P 

: mass density of air. 

CD 

CD 

wing torsion modes. 
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Commonly used superscripts 

1, 2, 3 : aerodynamic load parameters related to the 

various structural design requirements. 

' : three-dimensional airfoil aerodynamic quantity 

(i.e., a quantity corrected for finite span 


effects ) . 

Commonly used subscripts 

i as in A. : element i of vector A, or 
1 

: discretized value of function A(y), i.e., A^ = 

A(Yi). 

ij ; element ij of matrix. 

Special symbols 
{ } as in |A} : vector A, or 

: discrete values of function A(y) collected in 
vector |A]'^ = { A( y^ ) , A( y 2 ) , • • - } 

[ ] : square matrix. 


C.2 GLOBAL STRUCTURAL LEVEL WING MODEL 
C.2.1 Model Description 

To study the wing behavior at the global structural level, 
we model it as a straight beam aligned with the previously 
defined spanwise axis Y. We consider only two types of wing 


deformations . 
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i) Bending deformations in a direction normal to the 
wing mean aerodynamic chord, or, along axis N 
defined in Subsec. B.5.3. 

ii) Twisting deformations about the spanwise axis Y. 
This neglects chordwise bending. 

The parameters necessary to the description of the wing 
global structural behavior are the distributions of bending 
stiffness EI(y), torsional stiffness GJ(y) and total wing 
weight (structural and non-structural ) per unit length w(y). 
It must be noted that the wing is made of composite material 
as detailed in App. D. However the lay-ups are chosen so 
that no bending-torsion coupling exists. Therefore, the 
parameters retained will be adequate. 

The distributions will be written as follows 

El .P. (y) 

GJ.P. (y) 

w.P. (y) 

1 1 ' ' 


EI(y) 


GJ(y) 


■S 

n 

i=0 


w(y) 


n 




i=0 


166 



n (y-y^) 

j=0 ^ 



n 

s 


j=0 


n 

i=r> J 



i=0, . . . , n 


s 


(C.l) 


These expressions are interpolations based on Lagrangian 
polynomials (see Ref. 41). The EI^, GJ^ and are 

stiffnesses and weight values at selected, equally spaced 
stations y^ along the wing. They are the design variables 
of the global structural subproblem. They are collected in 
vector Xg and will be referred to as global stiffnesses. 

C.2.2 Elastic Axis and Center of Gravity Position 
The structural model described in the previous subsection is 
such that the elastic axis is assumed to lie on the spanwise 
axis Y. Clearly, the elastic axis position depends directly 
on the detailed wing design; it should be treated as a 
reverse control (see Chap. 3). However, this control will 
be ignored in a first approach to the design. It turns out 
that this is probably a very reasonable assumption, 
considering the detailed structural lay-out chosen (see App. 
D). 
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The position of a section center of gravity must be known 
to properly account for the inertia forces. This is also a 
reverse control. Initially, the mass center of gravity of a 
wing section will be assumed at the area center of gravity 
of that section, that is, at distances 6c from the section 
leading edge and nc above the reference chord (see Subsec. 

B. 2.2) . 

C. 2.3 Sailplane Weight 

The total sailplane weight is 


where includes the weight of the fuselage, tailplane, 
equipment and pilot; it is a fixed parameter in this design 
problem. For optimization and sensitivity analysis, the 


needed. If we perform the change of variables y=iib in Eq. 
C.l, we find 



derivatives of the weight with respect to b and the w^ are 


n 


s 


J=U 




n 



s 



therefore 
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C. 3 STRUCTURAL DESIGN REQUIREMENTS 

In order to obtain certification of an airplane by the 
appropriate authorities, one must insure that the structure 
is capable of withstanding a number of predetermined loading 
conditions. Also, static and dynamic instabilities should 
not occur within the flight envelope. Reference 40 
describes the requirements set by the Federal Aviation 
Agency (FAA) on sailplane designs. These requirements cover 
stresses in load cases corresponding to symmetric or non- 
symmetric flight conditions, towed flight and landing. Some 
requirements also cover flutter. In addition to satisfying 
the conditions set by the relevant regulating agency, the 
designer may have to investigate loading conditions specific 
to his own design. Clearly, a vast amount of calculation is 
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required in the structural design phase of a project. While 
modern optimization methods are perfectly capable of 
handling multiple design conditions, it is preferable to 
design the structure for a limited number of appropriately 
selected load cases and instability modes. The resulting 
design may then be checked for the remaining conditions and 
optimization may be performed again with additional load 
cases and instability modes, if desired. 

In this work, the structure will be designed to satisfy 
the three requirements described below. 

i) Requirement 1: limit out-of-plane bending. 

Sailplane wings are typically long and flexible. 

Excessive bending flexibility may result in control 

surface jamming, inappropriate ground clearance 

during landing or unsatisfactory dynamic behavior. 

We therefore require that the wing tip deflection 

be below a maximum limit u. at an indicated 

tmax 

airspeed = 85 m/s, and a wing load factor n^^^ = 
1. Since the wing does not present bending-torsion 
coupling, the constraint will control bending 
stiffness and weight distributions. 

ii) Requirement 2: limit stresses and preclude local 

instabilities. To size the details of the 


structure (see App. D), we require that no failure 
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occurs at a flight condition defined by n^^^ = 5.33, 
= 100. m/s. This corresponds to a pull-up at 
maximum authorized dive speed. This is a design 

condition, therefore; an ultimate safety factor 
= 1.5 must be used when calculating the ultimate 
loads . 

iii) Requirement 3: limit wing divergence dynamic 

pressure. The FAA specifies that the sailplane 
must be demonstrated "free of undesirable flight 
characteristics for speeds below the dive speed" 
(V^). The maximum authorized dive speed for this 
type of sailplane is about = 100 m/s. We 

require that the wing divergence speed Vj^ be such 
that Vj^ ^ 120. m/s. For straight wings with no 

bending- torsion coupling, divergence is torsional 
in nature. This constraint will therefore control 
the torsional stiffness distribution. 

C.4 WING TIP DEFLECTION CALCULATION 

C.4.1 Total Spanwise Distributions of Force and Moment 
At the end of the optimization of the wing aerodynamic 
design the spanwise distributions of aerodynamic force 
faN(y) moment m^Y(y) well as the total sailplane load 

factors n^j and n^ are calculated for the wing load factors 
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and flight speeds corresponding to Requirements 1 and 

2. Discrete values of force f and moment m . for 

aNi aYi 

various spanwise stations are transmitted to the global 

1 2 

structural subproblem in vectors A and A (see Subsec. 
B.5.4). The spanwise distributions of force and moment may 
be approximated by interpolation based on Lagrangian 
polynomials . 


n. 




i=0 

n.. 


^aNi^i(y) 


n>ly(y) = .E 


(y) 


Pi(y) 


^f 

n (y-y^) 

j=0 ^ 

xi 


y^ = ib/n^; i=0,...,n^ (C.4) 

To determine the total spanwise distribution of forces 


and moments, we add 

the 

contributions 

of 

inertia to 

the 

distributions given 

in 

Eq. C.4. 

The 

situation 

i s 


illustrated in Fig. 19 for a wing section at spanwise 
station y. The inertia forces are assumed to act at the 
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section area center of gravity. Measured positively along 
the axes C and N, the moment arm of these inertia forces 
with respect to axis Y are, for small wing twist (small 
E(y)-r) 


The distributions of transverse force and moment are, 
finally 


The distribution of chord c(y) is known to the global 
structural subproblem since it is communicated from the 
aerodynamic subproblem in vector H (see Fig. 4 and Subsec. 

B. 5.1) . 

C. 4.2 Wing Bending Response 

The wing bending model is a cantilever beam of varying 
bending stiffness EI(y) (Eq. C.l) under varying transverse 
load fjj(y) (E<3- C.6). In this subsection, we describe the 


e(,(y) = (6-6)c(y) 


ej^(y) = (n-y)c(y) 


(C.5) 





(C.6) 


calculation of the wing tip displacement using a single mode 
Rayleigh-Ri tz analysis. 
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^ Axis N 




Wing section, forces applied and corresponding 
geometrical parameters . 
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C.4.2.1 Wing Tip Deflection 
Assume that the displacement mode of the wing is given in 
terms of the non-dimensional spanwise coordinate Ti=y/b as 


u ( n ) = Du ( T) ) 


(C.7) 


where D is a constant. The specific expression chosen for 

"k 

u (n) will be given in the next subsection. By requiring 
that the work of the external force fjj(Ti) equals the bending 

■k 

strain energy for a virtual displacement about u (n), we 
find the tip displacement u^ as 


u^ = u(l) = u (l)b 


fjj('n)u (Ti)dTi 


EI(ti) [u*"(n) ]^dn 


(C.8) 


where ( )' = 9( )/3ti and ( )" = 9^( 


C.4.2.2 Assumed Displacement Mode 

The shape of the assumed mode will be taken as the 
displacement of a beam with linearly varying bending 
stiffness under linearly varying transverse load (hereafter 
referred to as model beam) . Using the superscript 1 to 
refer to the model beam, we have 

EI^(ti) = EI^(1 + &ti) 

f^(T)) = f^(l+an) 


(C.9) 
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The model beam obeys the following differential equation and 
boundary conditions 

( (1 + 3ti)u"(ti) ]" = kd+an) 


k 


1 4 

f b /El 
o 


1 

o 


u(n ) 


T1 = 0 


u’ (n) 1 


n=0 


0 


( 1 +N)u"(n) 1 ^^^ = I (1 + N)u"(ti) ] ' 1 ^^^ = 0 


(C.IO) 


The solution to Eq. C.IO is 
★ 

U ( T1 ) = kU ( T1 ) 

U*(n) = A^n'^+A3Ti^+A2Tl^+A^ii+Ag(l + &Ti)ln(l + BTi) 

A^ = -(a(-l+3&^+23^)+3&(l+e)^]/6B^ 

= [a(l-3&^)-3&(!+23) 1/123^ 

A 3 = (3&-a)/36P^ 

A^ = a/72B 

Ag = -A 3 /B (C.ll) 

★ 

The mode shape u (ti) in Eg. C.ll will be used together with 
Eq. C .8 to calculate the displacement of the wing tip. 

In order to carry out these calculations, we need to 
estimate the parameters a and & for the model beam. 
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Parameter a is found by requiring that the actual load fjj(Ti) 


A parallel approach was used to estimate parameter p. 
However, this occasionaly turned out to result in equivalent 
linear stiffness distributions that became negative. 
Therefore, parameter & is obtained by requiring that the 


have equal end values. This permits one to avoid negative 
stiffnesses, as the end stiffnesses EI(0) and EI(1) are 
design variables at the global structural level. They can 
be kept positive at all time during the optimization process 
by suitable side constraints. We have 


and the linear one f^(ri) are statically equivalent. That is 




This yields 


a = (3l2-1.5I^)/(I^-1.5l2) 



(C.12) 


actual distribution of stiffness EI(ti) and the linear EI^(ti) 
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This gives 


B = EI(1)/EI(0)-1 


(C.13) 


C.5 WING DIVERGENCE DYNAMIC PRESSURE 

A wing section is depicted in Fig. 20. If the wing was 
rigid, the angle of attack measured with respect to the 
section reference chord would be a^. The total twisting 
moment on the wing would be mY(y), as given by Eq. C.6. If 
the wing flexibility is taken into account, the section 
undergoes a twist 0 which results in an increment in the 
lift coefficient. 


This elastic lift causes an increment in torque per unit 
length (elastic torque) 


Cie(y) = c^^(y)0(y) 


(C.14) 


nieaY^Y^ = • SpV^c^ ( y ) [ 3-x ( y ) ] c^^ ( y )0 ( y ) 


or, if we set 


q = .5pvJ 


m 


eaY 


(y) = c^(y) iB-x(y) ]c^^(y) 


(C.15) 
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Figure 20: 


Lift coefficient induced by the wing elastic 
twist . 
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= qm*^Y(y)0(y) (C.16) 

* 

where q is the flight dynamic pressure. m , the elastic 

e a X 

torque per units of length, dynamic pressure and elastic 
deformation may be estimated on the basis of the information 

3 

contained in the vectors A {subvector A ) and H (see Fig 4). 


’"eay'y* 




m 


eaYi 


c^(yi)(e-xj 


)c 


' 3 
la i 


Pi(y) 


n (y-y. ) 

j=0 ^ 


f 


y^ = ib/n^; i=0,...,n^ (C.17) 

Vector H contains what is necessary to calculate c(y^) (see 

3 3 

Eqs . B.l and B.49), while vector A contains the x. and 
' 3 

c, ( see Eq. B. 66) . 
lai ^ ' 
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C.5.1 Wing Divergence Dynamic Pressure Calculation 
The torsional model of the wing is a cantilever beam of 
varying stiffness GJ(y) (Eq. C.l), under varying twisting 
moment m(y) = mY( y ) +m^^^( Y ) (Eqs. C.6 and C. 16-17). The 
divergence dynamic pressure can be calculated using a single 
mode Rayleigh-Ritz analysis. Using again the normalized 
coordinate t) , we write the wing twist as 

0(H) = E0*(ti) (C.18) 


[s 

quality 


For a given dynamic pressure, the wing twist amplitude E may 

be found by requiring that the work of the external forces 

m(Ti) equals the torsion strain energy for a virtual 

* 

displacement about 0 (n). 


b / m^ ( H ) 0 ( H ) dn 


E = 


GJ(H)0*’^(Ti)dTi - qb^ / 


(C.19) 


where, again, ( )' = 3( )/3 h- The divergence dynamic 

pressure is that which makes E unbounded 




GJ(T|)0*’^(Ti)dti 




(C.20) 


For this analysis, we will assume the torsion mode to be 
the divergence mode of a uniform wing with uniform twisting 
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moment per unit elastic twist (see Ref. 7). 


0 (ti) = sin(iTTi/2) 


(C.21) 


C.6 CONTROLS CALCULATION 


C.6.1 Wing Element Target Stiffnesses (Vector ) 


When performing the detailed structural design, the wing is 
modelled by n^ cylindrical wing elements of equal length and 
constant geometric and structural properties. Therefore, a 
wing element spans the distance between two consecutive 
stations used in the definition of the distributions of 
stiffness at the global level (see Subsec. C.2.1). Three 
global structural characteristics must be known for the 
synthesis of any one element; for element j, they are given 
by the vector 


These target flexural stiffness, torsional stiffness and 
weight per unit length will be taken for a given element as 
the averages between the corresponding global level 
variables . 



T 



(C.22) 


El . = .5(EI . t+EI . ) 
ei J-1 3 ' 


GJ . = .5(GJ. t+GJ.) 
ej j-1 3 ^ 


w 



(C.23 ) 
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The matrix is needed while developing the linear 

approximation to the constraint violations for local 
structural subproblem j. That matrix is easily obtained 
from Eq. C. 23 . 

C.6.2 Element End Forces (Vector F.) 

J jJ. 

The detailed structural design of each of the n^ wing 

elements is made to the values of the forces and moments 

calculated at the inboard side of the element. This seems 

justified in view of the type of loading chosen. 

The forces acting at the inboard section of element j are 

collected in vector F. 

1 

{F . = {V . ,M . ,T . 1 (C.24) 

I J J l 

with Vj , , Tj being respectively shear force, bending 

moment and twisting moment. To obtain { F^ | , the flight 

conditions specified for Requirement 2 (see Sec. C.3) are 

used to determine the corresponding spanwise distributions 

2 2 

of aerodynamic force moment (see 

Subsec. B.5.4); this is done in the aerodynamic subproblem. 
These quantities are passed in vector A to the global 
structural subproblem where they are modified to account for 
inertia effects (see Subsec. C.4.1). Finally, numerical 
integration is used to obtain the components of {F^}. 
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n 


V . 
1 


= I 


1=0 J^j 


- ^ w. P. 

^ i=o ^ 


(y)dy 


n. 


n 


M 




b , 


. , f „ . / P • ( y ) ydy 

2 ■ aNi /y . i w / j j 

1=0 J ^2 


4 1 r p. 

N 1 /y ■ 1 

1=0 J^2 


, (y)ydy 


n. 


= z 


b , 


m ^ . / P . ( y )dy 

1 . aYi Y . 1 ^ ^ 

^ 1=0 J^j 


n 


+ [n^(6-3)-n^(n-2r) ] 


^ “i jy^ '“i 


(y)c(y)dy 


y^ = (j-l)b/ng; j=o,...,rig 


(C.25) 


where y^ defines the inboard section of element j . 

The matrices 3F./3X , 9F./8A and 3F./3H are needed for 

J J 3 

the development of linear approximations to the measures of 
constraint violation in local subproblem j . The two former 


matrices can be easily obtained analytically from Eq. C.25, 
the latter will be calculated by finite differences. 


Appendix D 

LOCAL STRUCTURAL SUBPROBLEM 

The task of the local structural subproblem is to find the 
detailed design of the wing which corresponds to the 
spanwise distributions of weight and stiffnesses chosen in 
the global structural subproblem. The optimization is 
subjected to the additional constraints that the stresses at 
selected points in the structure do not exceed appropriate 
allowables and that the various wing components be free of 
buckling. 

The wing is divided into equal length cylindrical 
elements, each of which is optimized separately. The 
overall element dimensions (length, chord, thickness) are 
fixed by the aerodynamic subproblem. The global structural 
subproblem imposes target values for the constant weight per 
unit length and flexural and torsional stiffnesses. It also 
supplies values for the out-of-plane shear force, bending 
and twisting moments at the inboard section of the element 
for stresses and stability calculations. The design 
variables of the local structural subproblem are the 
detailed element dimensions. The stresses and the panel 
buckling loads are found using conventional design methods. 
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D . 1 NOMENCLATURE 


a 

A 

[A] 

b 

c 

d 

D 

E 

El 

e 

G 

GJ 


1 

M 


N 

[Q] 


s 


S 


: spar cap width. 

: cell cross-sectional area. 

: panel extensional stiffness matrix. 

: wing element length. 

: wing element chord. 

: rib spacing. 

: panel bending stiffness. 

: material extensional stiffness modulus. 

: wing element target bending stiffness. 

: wing element elastic axis position. 

: material shear stiffness modulus. 

: wing element target torsional stiffness. 

: vector containing overall dimensions for wing 

element j . 

: developed length of shell, or, panel dimension. 

: bending moment at element inboard section. 

: number of layers in composite laminates. 

; number of elements at the local structural level. 

: in-plane end load on element panel. 

: lamina stiffness matrix in panel axes. 

: coordinate along shell midplane. 

: composite allowable shear stress in material axes. 
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t : 

thickness of wing element, panels or composite 
layers . 

T : 

twisting moment at element inboard section. 

T’ : 

Eq . D . 16 . 

V : 

shear force at element inboard section. 

W, : 

target weight per unit length for element, weight 
per unit length for component i of wing element. 

X , X . : 

eg' X 

chordwise position of element center of gravity, 
same for element component i. 

X : 

material allowable strength along the fibers. 

y : 

spanwise coordinate. 

Y : 

material allowable strength across the fibers. 

z 

position of lamina midplane with respect to 
laminate midplane. 

Z , Z . : 

eg 1 

element center of gravity position above section 
chord, same for element component i. 

Z : 

shear strain. 

£ : 

normal strain. 

0 : 

layer orientation with respect to panel axes. 

V : 

material Poisson's ratio, or ratio between 
structural and non- structural weight for wing 
element . 

P : 

material weight density. 

0 : 

normal stress. 

I : 

shear stress. 
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wing section cell twist angle. 

Commonly 

Used Subscripts 

c : 

spar cap. 

or : 

critical (buckling) value of panel in-plane load. 

C, T : 

compressive or tensile allowable. 

f s : 

leading or trailing edge sandwich panel face 
sheets . 

fw : 

spar web sandwich panel face sheet. 

S : 

leading and trailing edge shell. 

V, T : 

shear flow induced by shear force or twisting 
moment . 

W : 

web sandwich panel. 

X, y : 

refer to panel axes. 

1-5, 6f, 

6r, 7f, 7r : end loads on various element 

components (Fig. 23). 

CM 

t — 1 

refer to material axes. 

a, 3 : 

refer to axis system defined for sandwich panel 
buckling analysis. 


Commonly Used Superscripts 


C : 

carbon material. 

e : 

equivalent . 

k : 

layer number in laminate. 

r : 

rib . 

S : 

fiberglass material. 

1 : 

transposed vector. 
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( )' : refer to quantities in material axes. 

D.2 LOCAL STRUCTURAL LEVEL WING MODEL 
D.2.1 Wing Construction 

The wing construction is typical of present day technology 
(see Refs. 42-43). It is a shell reinforced by a straight 
spar normal to the plane of symmetry of the sailplane (Fig. 
21). The spar caps are built of carbon fibers oriented 
spanwise and resist most of the bending moment. The shells 
as well as the spar webs are of symmetric sandwich made of 
fiberglass face sheets with foam core. The glass fibers are 
oriented at ± 45° with respect to the spanwise direction, 
thereby providing high shear stiffness. 

Control surfaces (flaps/airbrake and aileron) run the 
whole length of the trailing edge so that only the portion 
of the wing between the leading edge and the rear web, the 
torsion box, is capable of carrying load. Finally, 
fiberglass ribs are positioned chordwise to increase the 
buckling loads of the spar caps and the sandwich panels; the 
rib thicknesses will be held constant during optimization. 
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View A-A 


Figure 21: 


Wing structure layout. 
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D.2.2 Model Description 

To conduct the design at the local level, we divide the 

torsion box into n^ equal length, cylindrical elements (Fig. 

22). The chord length c,^ and the maximum section thickness 

t of each element are constant; together with the element 

length b they are communicated to the local structural 

subproblem from the aerodynamic subproblem in the vector 

(see Subsec. B.5.2). For each element, the design variables 

are the spar cap width a and thickness t^, the leading and 

trailing edge total thickness t and face sheet thickness 

s 

t^ , the front and rear spar web total thickness t and face 

X 5 V 

sheet thickness 't£^/ and the rib spacing d. Note that, 
technically, the thicknesses of the composite lay-ups vary 
in a discrete fashion since integer numbers of plies must be 
used. Similarly, the rib spacing is a discrete quantity. 
However, these will be treated as continuous variables to 
avoid difficulties while optimizing the design. 

D. 3 INTERNAL FORCE DISTRIBUTION 

A wing section under bending moment, shear force and 
twisting moment is depicted in Fig. 23. The orientation of 
M, V and T are consistent with the definitions in App. C. 
This section is concerned with the calculation of the 


In this appendix, we drop the subscript j for the 
quantities relating to element j (i.e., c^ ) . 


4 
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5 

Element 1 

2 

3 

4 




Leading and trailing 
edge shells 



Upper and lower 
spar caps 


Figure 22: Wing model for local structural design. 

Details of a wing element. Dimensions are in terms of chord 
c and thickness t (n^=5). 
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internal forces introduced in the wing section by that 
system of end forces. 

D.3.1 Bending Moment Effect 

The carbon fiber spar caps are about twenty times stiffer in 
the spanwise direction than the fiberglass shells (Ref. 42). 
It is therefore legitimate to assume that the entire bending 
moment is resisted by these spar caps. Furthermore, they 
are symmetric and their “thickness is small with respect to 
the distance that separates them. Therefore, we assume them 
to be subjected to a constant normal load per unit length 
approximately given by 

M 

= (D.l) 

a(t-2t. -t ) 
f s c ' 


D.3.2 Shear Flow Distribution 

Shear forces and twisting moments induce shear flows in the 
panels constituting the wing section. Simple equilibrium 
arguments (see Ref. 44) show that the shear flows are 
constant on panels not subjected to normal edge loads and 
linear on panels subjected to constant normal edge loads. 
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Figure 23: 


Internal forces induced in a wing section by a 
combination of shear force, bending moment and 
twisting moment. 
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Therefore, as summarized on Fig. 23, constant shear flows 
exist on the leading edge shell (^ 2 )/ the trailing edge 
shell (N^)/ the front spar web (N^) and the rear spar web 
(N^). The shear flows vary linearly between and on 

the upper spar cap and between and on the lower one. 

Four relationships may be written which express shear 
flow conservation at the junctions of the various panels. 


N^+N^-N 

N^+N^-N 

N^-N2+N 

N^-N2+N 

Clearly 


6r 


7r 


6f 

7f 


0 

0 

0 

0 


(T).2) 


= '^7r 


(D.3) 


Equilibrium of the upper spar cap in the spanwise direction 
yields 


N, -N,,- = a 
6r 6f 


dN, 


dy 


(D.4) 


Finally, the twist angle per unit length for any of the 
three cells is given by (see Ref. 44) 
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0 . = 


2A. 

1 


N( s )ds 


G( s)t( s) 
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(D.5) 


where the integral is carried around the relevant cell. 
is the cell cross-sectional area and s is a coordinate along 
the cell perimeter. G(s) is the material shear modulus and 
t(s) the wall thickness; both may be variable. We assume 
that only the sandwich face sheets carry shear. The 
equivalent shear modulus for the spar caps is obtained by 
specifying the compatibility of shear deformations for the 
various layers. 


t G*^+2t. G® 
^e _ __c f s 

t +2t^ 
c f s 


(D.6) 


Q 

G is the shear modulus for the spar cap carbon material. 
Since the carbon fibers are oriented spanwise. 



where 
G^ is the 
fibers of 
direction . 


(D.7) 

is the carbon shear modulus in material axes, 
shear modulus for the fiberglass material, the 
which run at ± 45° with respect to the spanwise 
Therefore (see for example Ref. 45) 


E 


s 

1 


2(l+v 


s 

12 


) 


(D.8) 
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where and fiberglass extensional 

stiffness modulus and Poisson's ratio in material axes. 
Using the subscripts L, T and S to designate respectively 
the leading edge, trailing edge and spar cell we have (see 
Fig. 24) 


= 


4G®A, 


■n^Il N^tn 


t,r t,. 

f S fw 


0T 


4G®A„ 


■N3IT N^t- 


t^ t^ 
f s fw 


2A, 


■<^6f"^6r>^ (N5-N4)t- 


S L' c 


(t^+2t^g)G 


2t. G“ 
fw J 


(D.9) 


The geometric quantities Aj^, A^, Ag, Ij^ and 1^ are 
calculated in Subsec. D.4.1. 

The relationships developed here will be used in the 
following subsections to determine the shear flows induced 
by a shear force and a twisting moment. 


D.3.3 Effect of Shear Force at the Elastic Axis 
A shear force at the elastic axis of a beam induces only 
bending deformations. Furthermore, the shear force is 
related to the bending moment by 


IS 
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Definition of parameters used in cell twist 
angle calculation . 
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dM 


V = - 


dy 

or, from Eq. D . 1 


dN, 


dy 


V 


a(t-2t^ -t ) 
' f s c ^ 


N 



OF POOR QUAUTY 


(D. 10) 


Equations D.2-4 and D.9-10 may now be used to set up a 
system of six equations in the six unknown shear flows 
induced by the shear force. Using the subscript V to 
emphasize the fact that the shear flows are induced by the 
shear force, we have 


^sv'^'sv-Nerv = ° 


^4V"^2V'^^6fV ° 


,,-N 


V 


6rV 6fV (t-2t^g-t^) 


t^ ^2V'^t^ ^4V ° 

f s fw 


-N.,„-r — = 0 


t^ 3V t. 5V 
f s fw 


(t +2t^ )G 
' c f s ' 


e^^6fV'^^6rV^'^_ s^^5V”^4V^ 

fw 


= 0 


(D.ll) 


The three last equations state the fact that the various 
cells do not rotate when V is applied through the elastic 


axis . 
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D.3.4 Elastic Axis Position 

Because we consider the shear force acting normal to the 
wing plane, we are only interested in the chordwise position 
of the elastic axis. It is found by stating that the sum of 
the torques induced about any point of the section by the 
shear flows of Eq. D.ll is exactly equivalent to the torque 
caused by the shear force itself (see Ref. 44). Therefore, 
taking the moments about point P in Fig. 25 a) 


D.3.5 Twisting Moment Effect 

The shear flows induced by a twisting moment may be obtained 
from the relationships developed in Subsec. D.3.2 with the 
added condition that the sum of the torques induced by the 
shear flows about any point of the section must now be 
equivalent to the applied twisting moment. By 
compatibility, the three cells must have the same twist 
angle. 



+N 

6fV 6rV 


) — ] 


2 


at 


(D.12) 




(D. 13) 


Also, there is no normal edge load on the spar caps. 


N 


1 


0 


N 


6rT 


N 


6fT 


(D.14) 
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N, 


a) Elastic axis position calculation, 
(only the shear flows contributing to the 
twisting moment about P are represented) 



Reference position for 
calculation in global 
structural subproblem 




force 


b) Corrected torque calculation. 


Figure 25: Elastic axis position and corrected torque 

calculation . 
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where the subscript T indicates twisting moment induced 
shear flows. This yields the following system of five 
equations in five unknowns. 


^3T'^^5T"^6fT ° 


^4T"^2T'*'^6fT ° 




N 


t. 2T t. 4T A„t. 
fs fw T fs 




4aA^G‘ 


tA„ 


*^3T 


■^5T~ 


-N. 


'fs 


'fw 




(^5T ^4T^ 


= 0 


where the last equation states the equivalence of the 
torques. It may be deduced from Eqs . D.12 and D.14, 

replacing eV by T. 


D.3.6 Combined Loading 

The applied values of the shear force, bending and twisting 
moment to use in the design of a given wing element are 
supplied by the global structural subproblem. However, the 
twisting moment is calculated about a point located at the 
maximum airfoil thickness, halfway between the upper and 
lower spar surfaces. The twisting moment must be corrected 
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to account for the offset of that point with respect to the 
elastic axis. To obtain the internal forces corresponding 
to any set of applied forces, the following steps must be 
carried out. 

i) Obtain the normal edge load on the spar caps from 
Eg. D. 1 . 

ii) Calculate the shear flows induced by a shear force 

assumed to act at the elastic axis using Eq. D.ll. 
The now known. 

iii) Find the elastic axis position from Eq. D12 . 

iv) Calculate the total moment of the applied forces 
about the elastic axis (from Fig. 25 b) 

T’ = T+V(|-e) (D.16) 

v) Calculate the shear flows induced by the complete 
twisting moment using Eq. D.15, and obtain the 

vi ) Combine the N^y from step ii) with the from 

step v) . 

D. 4 WING SECTION GEOMETRIC PROPERTIES 

D.4.1 Shell Developed Lengths and Cell Cross-Sections 
Referring to Figs. 22 and 24, we have 

= .25tc-.34ta 


ta 
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= ,27tc-.35ta (D.17) 

= [ ( ,3t+.02c)^+( ,37c-.5a)^]^/'^+.35t 

+[ ( .35t-.02c)^+( .37c-.5a)^]^/^ 

= 2[ ( .3t)^+( ,38c-.5a)^]^/^+.4t (D.18) 

D.4.2 Element Weight Data 

The calculations for the element weight per unit length and 
the position of its center of gravity are detailed in Table 
2. The position of the various items are given in an X-Z 
system of axes centered at the chordwise position of maximum 
airfoil thickness, halfway between the upper and lower 
surfaces (see Fig. 22). The item numbers are defined in 
Fig. 26. 

The entries in Table 2 are in terms of the variables 

SC f 

defined in Fig. 22. p , p and p are the weight densities 
of, respectively, the fiberglass sandwich face sheets, the 
carbon fiber spar caps and the sandwich foam core. Item 11 
represents the total weight per unit length of the ribs. 
The wing element is assumed to have b/d ribs, the weight of 
which is proportional to the torsion box cross-sectional 
area and the weight per unit area p^. 


Finally, item 12 is 
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TABLE 2 

Weight and center of gravity position for wing element 

components . 


Item 

w . 

z . 

X . 


1 

1 

1 

1 

[ ( .37c-.5a)^+( .3t+.02c)^]^/^T* 

.35t-.01c 

- . 185c- . 25a 

2 

.35tT 

. 025t- . 02c 

-.37c 

3 

[ ( .37c-.5a)^+( .35t-.02c)^]^/^T 

- . 325t- . 01c 

- . 185c- .25a 

4 


.0 

- . 5a 

5 


.5t 

.0 

6 

a(t^P*^+2t^gP®) 

-.5t 

.0 

7 


.0 

. 5a 

8 

[ ( .3t)^+( .38c-.5a)^)^/^T 

. 15t 

. 19c+ .25a 

9 

[ ( .3t)^+( .38c-.5a)^ 

-.15t 

. 19c+ . 25a 

10 

■4tl2t^^p=MV2tj^)p^l 

.0 

. 38c 

11 

(Aj^+Ag+A^)pVd 

.0 

.0 


11 



12 

V Z-. w. 

.01c 

.46c 


i=l ^ 
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11 Wing ribs 

12 Non- structural components 


Figure 26: 


Numbering scheme for element weight and center 
of gravity calculations. 
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the total non- structural weight per unit length (control 
surfaces, control linkages, . . . ) . It is assumed to be a 
fixed percentage of the structural weight (v) and to be 
concentrated far behind the spar. The total element weight 
per unit length is 


12 




i=l 


Its center of gravity is defined by 


(D.19) 


X 


eg 


z 

eg 



x.w. 
1 1 


z . w . 
1 1 


(D.20) 


D.4.3 Element Bending Stiffness 

Assuming that the entire bending stiffness is provided by 
the spar caps we have, approximately 


E^^at E®at 

El = r-^(t-2t. -t — ^-^[(t-t. )"^+(t-3t. -2t y] 

2 ' fsc' 2 ' fs' ' fs c' 


(D.21) 


where the last term is the contribution of the fiberglass on 
either side of the spar caps. E , the carbon fiber spanwise 
extensional stiffness modulus, is 
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(D. 22 ) 


where is the material principal extensional stiffness. 
Since the glass fibers are oriented ± 45° with respect to 
the spanwise direction. 


E^ = 4 




12 


(D.23) 


where E^, ^^2 fiberglass principal 

extensional stiffness, shear stiffness and Poisson's ratio, 
in material axes. 


D.4.4 Torsional Stiffness 

The torsional stiffness of a beam is defined as the ratio 
between an arbitrary applied torque and the resulting twist 
angle. By compatibility, the twist angle of a multicell 
beam is the same as that of any of its cells. 


GJ = 


2A.T 

1 




N( s)ds 


/^G(s)t(s) 
or, for the leading edge cell 


GJ = 


2A^G®T 

Li 


^4T^^ ^2T^L 


2t- 2t. 

fw f s 


(D.24) 


(D.25) 
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The torsional stiffness is obtained by solving the torsion 
problem (Eq. D.15) with an arbitrary torque and then using 
Eq . D . 25 . 


D.5 CALCULATION OF STRESSES 
D.5.1 General Procedure 

The different fiber reinforced laminated panels constituting 
the wing torsion box are of symmetric lay-ups. Furthermore, 
they are assumed to be loaded by in-plane edge loads. 
Therefore, their deformations entail only in-plane strains. 
For each laminated panel (laminate) the strains are related 
to the edge loads by® 


\t] = (A]‘^|N1 


[t] = { E , E , y j 

{ j * X y xy’ 


{N} = {N ,N ,N } 

j j * X y xy‘ 



n. 


z 

k=l 




i, 3=1,2, 6 


(D.26) 


The strain and load vectors (e,N) are given in a set of axes 
called laminate axes (see Fig. 27). The laminate 
extensional stiffness matrix A is obtained from the 


5 


The developments of this section belong to the classical 
lamination theory, they are taken from Ref. 45. 
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transformed reduced stiffness matrix Q and the thickness t 
of the n^^ layers constituting the laminate. The transformed 
reduced stiffness matrix of a layer is dependent on the 
characteristics of the material making up the layer and on 
the orientation of the layer material axes with respect to 
the laminate axes. General expressions for the Q matrix 
may be found in Ref. 45. 

The stresses in a layer are given in laminate axes by 




. k,T 

.0 I = 


, k k k , 

to ,0 , T } 

‘ X y xy’ 


(D.27) 


Finally, the stresses in the layer material axes are found 
from 


{o'^l = [T^] jo^] 



II 

Q 

k k 
°2''^12 

1 

IT) = 

r k 2 

c 

k2 

s 

k k 
2s c 


k2 

s 

k2 

c 

k k 
-2s c 


k k 

hs C 

k k 
s c 

k2 k2 
c -s 


k ok k ok 

c = cos0 ; s = sin0 


(D.28) 


where the ' symbol is used to identify quantities measured 
in material axes. The stress analysis for the panels will 
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Laminate axes: XY 

Material axes for layer k: 12 


Figure 27: 


Definition of axes system and applied edge loads 
for an arbitrary laminated panel. 
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be discussed in the next two subsections. The necessary 
information is summarized in Table 3. 

D.5.2 Stresses in the Spar Caps 

The details of the spar cap construction are depicted in 
Fig. 22, while Table 3 gives the location of the points to 
be monitored (points 1 through 8). The transformed reduced 
stiffnesses are now given for the various layers (see Ref. 
45). For the carbon fiber material 













(D.29) 



material extensional stiffness moduli in the fiber 
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TABLE 3 

Summary of information for stress analysis. 
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and across the fiber directions, Poisson's ratio and shear 
stiffness modulus. For the fiberglass material (thickness 


2t. ) 

f s' 


^11 = ^22 


+G 


2(1-v“2) 


•12 


=66 


2(1-v,2) 


2(1^v-2) 


-G 


12 


12 


2i6 = 226 = ° 


(D.30) 


where the superscript s refers to the fiberglass material. 
The stresses in the various layers may be determined from 
Eqs. D. 26-27, together with Eqs. D. 29-30. They are in 
material axes for the carbon fiber layers; for the 
fiberglass layers, however, transformation D.28 must still 
be used. 


D.5.3 Stresses in the Sandwich Panels 

A uniform state of stress exists in the sandwich panels 
constituting the leading and trailing edge shells, as well 
as the spar webs (see Fig. 23). Hence, the stresses need 
only be monitored at the four points mentioned in Table 3 
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(points 9 through 12). The sandwich core will be assumed 
non-structural . In this case, the calculations of Subsec. 
D.5.1 may be carried out quite simply. For the leading and 
trailing edge shells, we have, in material axes 
N 

_ _ xy 

°1 ~°2 2 t . 

f s 

1^2 = 0 (D.31) 

while for the spar webs, t^^ in Eq. D.31 should be replaced 


D.5.4 Failure Criterion 


The states of stress described in the two previous 


Therefore , 


failure 


subsections are two-dimensional, 
criterion must be used to determine whether a given loading 
induces failure. Tsai-Wu's criterion reads (Ref. 45)® 

.2 

’ (D.32) 


Oi(Oi.Xc-Xt) 02(02 -Y^-Y,^) 




•+— -1 < 0 


Vc 


Y Y 
T C 




where and X^ are the material compressive and tensile 

strengths along the fibers, Y^ and Y,^, the corresponding 
quantities across the fibers and S the in-plane shear 


® General formulations for Tsai-Wu's criterion usually 
involve a strength biaxial tension test. Such a stress is 
difficult to determine accurately. In the formulation 
chosen here, the need for that strength is alleviated by 
specifying that Tsai-Wu's criterion reduces to Von Mise's 
yield criterion when applied to an isotropic material. 
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strength. Failure will not occur in the wing section if 
inequality D.32 is satisfied at all points defined in Table 
3; there are twelve such conditions. 

D. 6 SECTION STABILITY 

The stability of the wing box will be ascertained by 
verifying the stability of each individual component. The 
leading edge shell, which exhibits a strong curvature, will 
be replaced by two flat panels. because of the model cross- 
section symmetry and the nature of the system of induced 
loads (Fig. 23), only seven panels need be considered; they 
are numbered as shown on Fig. 28. Two different buckling 
problems must be solved: 

i) Buckling of orthotropic sandwich panels under in- 
plane shear (panels 1 through 5). 

ii) Buckling of orthotropic composite panels under in- 
plane compression and shear (spar caps, panels 6 
and 7 ) . 

These problems are discussed in the following subsections. 
D.6.1 Sandwich Panels Buckling 

The sandwich panels constituting the wing torsion box are 
primarily loaded in shear. This subsection is concerned 
with the s calculation of the critical load of orthotropic 
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Figure 28: 


Panel numbering for cross-section stability 
analysis . 
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sandwich panels; the developments are adapted from Ref. 46. 

Figure 29 defines the nomenclature used. Note that the 

panel edge of length 1 is that which is shorter. The glass 

a 

fibers are oriented at ± 45° with respect to axes o and &. 
In this set of axes, the face sheets material properties are 



4 

[2(1-v®2)/E®+1/G®2] 



E 


s 

1 


2 ( 1+v 


s 

12 


) 


U/G®2-2(1-v®2)/E®] 

= (D.33) 

[1/Gi2"2(1-v^2)/Ei1 

s s s 

where E^, G ^2 '^12 fiberglass material 

properties in material axes. Neglecting the contributions 
of the core material, the panel bending and twisting 
stiffnesses are 


= 


E^[t^-(t-2t^)^ 






aB 


12 


(D.34) 
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Figure 29: Nomenclature used in sandwich panel buckling 

analysis . 
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N 

The critical shear flow is 


4jD 

a 

Kj^ and Kj^q are given by 


^MO 

A 


C^+2B2+B2+AV( 1+1/C^) 
1+7(03+83+ (0^+83 )/C^]+V^A/C^ 

= Kj^(V=0) 

= 1-B^*B3(Cj+2B2+C3) 


^2 = 

®3 = 

Cj = = I/C3 = 

For an isotropic core material 


l^(t-2t.)G^ 
a I 


(D.36) 


(D.37) 


being the core material transverse shear stiffness. 
Finally, j is a function of and 1^/1^; for simply 

supported panels, it is approximately given by 
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+(3. 1502-8. 5488B„)(1 /I.) 

2 a p ' 


3 


+ (-. 4501 + 5. 1740B2)(l/l3) 


4 


(D.38) 


The critical load (Eq. D.35) must be calculated for the five 
different sandwich panels. The necessary information is 
summarized in Table 4. For stability, we require for each 
panel 


xycr 

D.6.2 Spar Caps Buckling 

The spar caps are loaded in their plane by spanwise normal 
forces and shear forces. We now calculate their critical 
loads following developments given Ref. 47. We make the 
conservative assumption that they are very high aspect ratio 
simply supported plates. We assume further that the shear 
force is constant along the edges, taking as applied shear 
the largest of and (Sec. D.3). The notation used is 

that of Sec. D.5. 

The spar cap bending stiffnesses are 


N 


N 



(D.39) 


n 


1 



(D.40) 
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Summary of information for sandwich panel buckling analysis. 


Panel 

N 

xy 

1 1 ^ ^ 
a' e 

t 

^f 

1 

^2 

d, [ ( .37c-.5a)^+( ,5t)^]^/^ 

t 

s 

^fs 

2 

N_ 

d, .4t 

t 



3 


s 

fs 

3 

^3 

d, [ ( .38c-.5a)^+( ,3t)^]^/^ 

t 

s 

f s 

4 

Nc 

d, t 

t 



5 


w 

fw 

5 

N/i 

d, t 

t 



4 


w 

fw 


( * ) 

1^ should be taken as the smallest of the two values 
given, 1^ is the remaining one. 
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The Q? . are the transformed reduced stiffnesses (Subsec. 

k k 

D.5.2) for layer k, t is that layer thickness, and z the 

distance between that layer midplane and the laminate 

midplane . 

The critical compressive load is given as 


N 

xcr 






1 + 


(^11^22) 


1/2 


(D.41) 


where c is the length of the edge under the normal loads. 

The critical load in shear is given in terms of the 
parameter 


^^iA2 \ 
»? 12 '^^^ 66 / 


(D. 42 ) 


Then, for 0 S 6 < 1 


N 

xycr 


4-||D22(Di2-2D^g 

c 


1/2 


C = 11. 71+. 0950+1. 7670^ 
a 


(D.43) 
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For panel stability under combined loading, we require 



(D.45) 


This condition must be satisfied for both spar caps, the 
information necessary to conduct the calculations is 
summarized in Table 5. 
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TABLE 5 

Summary of information for spar cap buckling analysi 


Spar Cap 

N 

X 

N 

xy 

c 

upper 

-^1 

max(Ngf,Ng^) 

a 

lower 

^1 

max 

a 
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NUMERICAL DATA 

This appendix gives numerical values for the parameters 
introduced in Apps. A-D. These values were used when 
solving the two- level structural design problem or are 
appropriate for completion of the four-level problem. 

E. 1 PERFORMANCE SUBPROBLEM 

The cross-country speed may be calculated for any of the 
three thermal models described in Fig. 13; the rate of climb 
constraint should require a minimum rate of climb value 
~ • ”7 m/sec, in a weak thermal. An interthermal 

cmin 

downdraft of = .3 m/sec appears typical. The mass 

density of air is p = 1.225 kg/m^ (standard day, sea level), 
the acceleration of gravity is g = 9.81 m/sec^ . 

E.2 AERODYNAMIC SUBPROBLEM 

The kinematic viscosity of air is v = 1.4607 10“ ® m^/sec 

(standard day, sea level). The wing angle of attack is 
taken as i^ = .0349 rad (2°). The following sailplane 
parameters were mostly estimated from sketches of the Nimbus 
II sailplane in the references to App. A. The sailplane 
fuselage is characterized by a radius b^ = 0.4 m, a maximum 
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N 

height h^ = 0.8 m, a total length 1^ = 7.5 m, a maximum 
cross-sectional area A = 0.6 m^ , a wetted area . = 10.0 
m^, a planform area = 3.0 m^, a distance nose-wing 1^^ = 

2.0 m, a total volume = 1.3 m^, a distance between wing 
and horizontal tailplane 1^^ = 5.0 m, the wing planform area 
inside the fuselage is = 0.6 m^ . Reference 39 (Chap. 
IX) recommends a static stability margin (x^-x^^)/c = .05. 
The horizontal tailplane has a planform area = 1.2 m^, an 
average chord = 0.3 m, an aspect ratio = 5.0, a taper 
ratio = 0.5, a thickness ratio (t/c)j^ =0.1, a vertical 
distance between wing and tailplane hj^ = 1.2 m. The 
vertical tailplane has a planform area = 2.4 m^ , an 
average chord c^ = 0.5 m, a thickness ratio (t/c)^ = 0.1. 
The two-level subproblem is solved using a fixed total 
sailplane weight W = 4000. N for the aerodynamic 
calculations; the wing is assumed untwisted = 0), 
the planform is described by b = 10.0 m, b^^ = 5.0 m, c^ = 

1.0 m, c^ = 0.7 m, and c^ = 0.3 m. The Fourier series 
analysis for spanwise distributions of lift and drag is 
based on n^ = 5 control points; while the aerodynamic 
parameters transmitted to the global structural subproblem 
(Subsec. B.5.4) involves 6 control points (n^ = 5). 
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E. 3 GLOBAL STRUCTURAL SUBPROBLEM 

In addition to the parameters related to the structural 
design Requirements and given in Sec. C.3, it is specified 
that the tip displacement never exceed 0.7 m. The weight 
for the sailplane without wings is W = 2000. N. The 

discretization of the distributions of weight per unit 
length and bending and torsional stiffnesses are based on n^ 
= 5 intervals. The following side constraints are used 1.0 
< Wj < 1.0 101“ (N/m), 1.0 10^ < Elj < 1.0 10 1 “ (Nm^), 1.0 

1Q3 < GJ . < 1.0 lOi " (Nm2 ) . 

E. 4 LOCAL STRUCTURAL SUBPROBLEM 

The discretization at the lowest structural level involves 

n , = 5 elements, 
el 

The unidirectional carbon fiber material is a high- 
Strength graphite-epoxy with the following parameters : 

= 1.47 1011 N/m2, = 1.09 10i“ N/m^ , 

N/m^, v^2 = 0.38, = 1.41 10^ N/m^ , = 1.46 10^ 

N/m^, = 1.48 10® N/m^ , = 4.2 lO^ N/m^ , S = 9.52 

10'' N/m^, = 1.58 10^ N/m^. 

The fiberglass material is a bidirectional cloth with 
fibers running in perpendicular directions and having the 
following properties: = E^ = 1.55 lOi ° N/m^ , G®^ = 

5.5 10® N/m2, v^2 = -25, = 8.0 lO^ N/m^ , xj = yJ = 

Xj = Y® = 2.4 10® N/m2, p® = 1.67 10® N/m^. 
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The foam used as core for the sandwich panels is 
described by = 4.0 10® N/m^ , = 5.89 10^ N/m^. The 

wing ribs are constructed of a 5.0 mm thick foam core with 
two fiberglass face sheets of 1.0 mm each, for a total 
weight per unit area p =36.3 N/m^ . The ratio between non- 
structural and structural wing element weight is v = 0.5. 

The design variables are subjected to the following side 
constraints: (all dimensions are in meters) 5.0 10'^ ^ a < 

1.0 10'^, 2.0 lO"* 5 t < 5.0 10-2, 2.5 lO'^ ^ t < 5.0 

o s 

10-2, 5 0 10-^ < t. ^ 1.0 10*2, 2.5 10-2 ^ t S 5.0 10'2, 

fs w 

5.0 lO-** ^ t. ^ 1.0 10-2, 5 0 10-2 ^ d ^ 5.0 10-1. 

fw 

Furthermore, the following constraints are added to preserve 
the validity of the analyses: a < 0.1c, t^^ S 0.2t^, t^^ < 
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